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ABSTRACT 


We address critical issues arising in the practical implementation of processing real 
point cloud data that exhibits irregularities. We develop an adaptive algorithm based 
on Learning Theory for processing point clouds from a stationary sensor that stan- 
dard algorithms have difficulty approximating. Moreover, we build the theory of 
distribution-dependent subdivision schemes targeted at representing curves and sur- 
faces with gaps in the data. The algorithms analyze aggregate quantities of the point 
cloud over subdomains and predict these quantities at the finer level from the ones 


at the coarser level. 
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INTRODUCTION 


Point clouds are sets of unstructured points with no specific ordering or connections. 
In typical settings relevant to several applications, the data points correspond to a 
surface. For instance, terrain data has that the data points represent a location on 
earth. The data points are received by processing images, video, manual measure- 
ments, etc. One of the most promising ways of receiving terrain data is through 
LiDAR, which is a remote sensing technology that measures distance by illuminat- 
ing a scene with a laser and analyzing the reflected light. Research problems often 
associated with this type of point clouds include: efficient feature extraction, surface 
reconstruction, and visualization. Several approaches such as subdivision, implicit 
surfaces, and variational surfaces have been developed toward these goals for point 
clouds with nearly uniform sampling at a rate much higher than the required reso- 
lution for reconstruction. When subjected to a relatively low sampling rate, LiDAR 
scanners commonly produce point clouds with undesirable properties such as miss- 
ing data, occlusions, or uneven distribution. To address these issues, we introduce 
mathematical theory and tools for processing irregularly distributed data. 

The content of this dissertation addresses two main problems. The first is the 
processing of point clouds that are regularly distributed in a topology different than 
the one provided by the standard Cartesian grid. The second problem is the approx- 
imation of distribution-dependent quantities during subdivision. A new paradigm is 
introduced for building subdivision schemes that take into account the distribution 
of the measure, thereby yielding schemes tailored to unevenly distributed data. 


The first problem is focused on surface reconstruction and starts with learning 


in a topology different from the Cartesian space. We assume the point clouds are 
sensed from a surface and are functional in the topology of the scanning process. In 
particular, we consider spherical learning, which is a natural approach for analyzing 
fixed or slow moving LiDAR scans. Fixed scans collect data in terms of two angles 
and the distance from the sensor which result in nearly regular distribution of points 
in spherical coordinates, yet the objects being sensed are mostly rectangular. In 
Cartesian coordinates, these point clouds are very irregular, as some regions are 
dense and others are sparse. Our algorithm partitions the surface in spherical space 
then builds piecewise linear approximations over the partition in Cartesian space. 
The first step is to transform the point cloud to the underlying topological space (in 
this case spherical) and generate a family of partitions using newest vertex bisection, 
which is organized in a multiresolution tree structure. The error is measured in the 
standard topology (in this case Cartesian). Thus we fit a Cartesian plane on each 
cell of the partition using the aggregate quantities of the cell. It is important to note 
that the basis of the algorithm calculations is the aggregate quantities of the point 
cloud data and not the individual data points. The aggregate quantities ought to 
characterize the local behavior of the data. We consider the moments of the point 
cloud within a subdomain as the aggregate quantities because the moments are useful 
tools in understanding the behavior of the points. Since we use the quantities and 
not the individual data points, we can easily process streaming data. When new 
points are added to the point cloud, we adjust the numerical values of the aggregate 
quantities and do not need to introduce a new object. Moreover, we develop a new 
method for encoding the surface approximation by introducing nonlinearity into the 
multiresolution tree itself. 

Our algorithms are based on the theory and methodology in Learning Theory 
by Binev, Cohen, Dahmen, DeVore, and Temlyakov (see [4] and [2]). These papers 


focus on the approximation of the surface constructed over the partition cells. Their 


work has provided the theoretical framework and tools from which we build our 
algorithm. Some explanation of the Learning Theory estimates in [4] and [2] as well 
as a description of the results will be given, but we will not repeat the estimates as 
they are very similar to what has been done and there is no added contribution for 
our particular properties. The estimates for piecewise constant approximations are 
optimal. The authors achieve similar results for piecewise polynomial approximations 
when a regularity condition is imposed on the distribution. Our practical algorithms 
use techniques such as singular value decomposition (SVD) to ensure that the process 
is stable. 

The second problem of distribution-dependent subdivision schemes focuses on the 
issue of visualization. The film and video game industries, for example, often employ 
subdivision as a prominent technique for efficient visual representation of comput- 
erized objects. Subdivision surfaces connect the ideas of continuous and discrete 
models. Notable examples of subdivision schemes include: Catmull - Clark [7], Doo 
- Sabin [9], and Loop [12]. In Section 1.4 we discuss the one-dimensional example of 
subdivision based on splines and discuss the two-dimensional example of 4-8 Subdi- 
vision by Velho and Zorin [15]. Most of this work concerns subdivision schemes for 
dense and evenly distributed data and is driven by the functional values. They have 
a very different flavor than the methods we discuss. Our approach in Chapter 3 is 
to approximate not only the surface but also the distribution-dependent quantities, 
which we define to be the desired quantities that depend only on the location and 
distribution of the points in the data set. As in the first problem, the moments are 
the local aggregate quantities and drive the development of the subdivision schemes. 
We draw analogies from polynomial reproduction to reproduction of the aggregate 
quantities of particular types of measures. This method of learning the distribution- 
dependent quantities is a novel method, so we start by building the foundation. This 


dissertation presents algorithms and results for curve subdivision over intervals and 


surface subdivision over triangles. 

The dissertation is organized in the following manner. In Chapter 1, we supply the 
definitions and background information used in the remaining chapters. In Chapter 2, 
we tailor the general algorithms in Learning Theory [4] and [2] to the case of spherical 
underlying topology. A nonlinear procedure to encode our surface approximation is 
described. Chapter 3 explores the idea of learning the measure to adjust for gaps and 
unevenly distributed point clouds. The methods explained in Chapter 3 serve as a 
launchpad for future work. Finally, Chapter 4 summarizes our results and outlines 
some of the many possible directions for future study. 

This work was partially supported by NSF grants DMS-0915104 and DMS-1222390, 
MURI/ARO grant W911NF-07-1-0185, and DEPSCoR ONR#N00014-07-10978. 


CHAPTER 1 
PRELIMINARIES 
1.1 POINT CLOUDS 


A point cloud D is a finite set of d+ 1 dimensional data points. Each point p € D 


is of the form p = (x,y) where x € R¢ and y € R. We assume that there are 


bounded regions X C R¢ and Y C R such that for all possible point clouds D and 


points p € D we have that x € X and y € Y. For convenience we assume X is a 


d-dimensional cube, as for every bounded region there is a hypercube containing it. 


We purposefully do not consider a point p to be a single vector in R“! because we 


assume that p is the result of probing a real valued functional surface, namely s, at 
x. Sometimes we will refer to y as the value at x, and at times it is representative of 
the height of a terrain at the location x. In the case of perfect acquisition, we would 
have y = s(x), but in reality we expect to have errors in acquiring both the value y 
and its position x. 

Mathematically, the point cloud is much more than a simple set of points. Here 
we establish the mathematical framework for the distribution of the points and the 
relation between the point cloud and the surface. We assume that D consists of in- 


dependent random observations that are identically distributed according to a prob- 


ability measure p on X x Y C R*!. Probability measures are the same as the more 


general notion of positive measure with the added condition that the measure of the 


entire probability space must be equal to one. The marginal probability measure on 


X is denoted by px and is defined for S C X as 
px(S) := p(S x Y). 


When extending X to a d-dimensional cube, we also extend px with zeros outside of 
the original set. The measure px is considered unknown. In Chapter 2 we assume 
that the underlying topology of px is spherical, and in Chapter 3 we approximate 
px on subdomains each time we subdivide. We let L2(X, px) denote the space that 
consists of all functions from X to Y which are square integrable with respect to 
px. We also have that dp(x,y) = dp(y|x)dpx(x), where p(y|x) is the conditional 
probability measure on Y given x € X. The regression function is the conditional 


expectation of y at x 
p(x) := | ydelyle). 
In this dissertation, we assume that there is an M such that |y| < M almost surely, 


which implies that |s,| <M almost everywhere with respect to py. Furthermore, s, 


is the function in D2(X, px) that minimizes the risk functional 


Joy ly FOOPde = lly — FO) exexex) 


(see, [3]). In this model we consider p to be a measure concentrated around the 
surface of interest, s, and will assume that s, is the surface that represents the point 
cloud. Ideally there is no measurement error in the acquisition of the point cloud, 
so Ss, = 8s. Therefore, we disregard the model error ||sp — s||r,(x,9,) in our future 
considerations and our objective is to approximate s,. We denote our approximation 
of s, by &. 

The points in D will be partitioned into subdomains. Throughout the paper we 
calculate important aggregate values from the group of points in each given subdo- 
main. The main idea is that on a subdomain R, we only examine these aggregate 


quantities representing the points in R and ignore the individual points. Let Q(R) be 


the aggregate representation of the points in the subdomain R. These aggregate quan- 
tities should be quantitative measures of the shape of the points in the subdomain. 


The moments of px are a natural choice to include in our aggregate representation. 


When d = 1 and x € R the moments are: 


M,(R) = a px ye fore 0152 Bw (1.1) 


xER 


We will also consider the moment M,(R) = i. y dpx. From these moments several 
useful quantities can be calculated, namely the average values and variance. The zero 
moment gives us the size of the measure px on R. M,/Mp gives us the average value 
of x on R, and M,/Mo is the average value of y on R. The variance of x on R is 
M,\? 
mcrae 
The moments for x = (x1, 22) € R? are: 


M,j(R) = [ai dpx © Se ai xh fora 0012 She (1.2) 


xER 


Expanding on the values discussed for R, we have the covariance which is Mq, — 


Mi oMo1. With these quantities we can build the covariance matrix. The smallest 
eigenvalue of the covariance matrix is another valuable quantity to be used in some 
algorithms. 

A quantity q is additive if q(R,) + q(R2) = ¢(Ri U R2) for disjoint subdomains R, 
and Rj. We desire that all of the quantities in Q are additive. If a quantity in Q is 
not additive, then future calculations will require more computations. For additive 
quantities, knowing the quantities over finer subdomains allows for quick calculation 
of the coarser subdomain quantities. For our purposes, Q(R) consists of the moments, 
and from the moments we calculate the other quantities as necessary. We choose the 
moments because they are additive and descriptive of the points they represent. We 
do not let the Q(R) be the averages, variance, etc. because these quantities are not 


additive. 


Point Clouds from LiDAR 


The point clouds considered in Chapter 2 are LiDAR scans, where the sensor remained 
in the same or nearly the same position. These point clouds generally are irregularly 
distributed in Cartesian space. For instance, the regions where objects are sensed 
are relatively dense and the regions without a sensed object are relatively sparse. In 
spherical space, these types of point clouds generally are evenly distributed with the 
exception of gaps or holes in the data. 

Two common issues from real world LiDAR scans cause gaps in the point cloud. 
The first is occlusion, which occurs when one object blocks the sensor from acquiring 
points in a region. The point cloud will only contain data along the line of sight of the 
sensors. In urban areas, it is common to have only points over one or two sides of a 
building, the other sides being occluded. The second issue is no returns. Sometimes, 
no information is transmitted back to the sensor, and thus no return is recorded. 
Where no returns occur, holes will appear in the point cloud. One potential cause of 
no returns is water. Infrared beams used by LiDAR tend to be absorbed by water 
and show very weak or no return signal when they hit rain, snow, thick smoke or very 


thick haze. 


1.2. BARYCENTRIC COORDINATES AND BERNSTEIN POLYNOMIALS 


A simplex is the convex hull of d+ 1 vertices v; € R?, i = 0,...,n, and is denoted 
by A = (vo, Vi,-.-,Va), a (d + 1)-tuple of d-dimensional vectors. If {v; — vo}“, 
are linearly dependent, then the simplex is degenerate. We want to emphasize that 
points will only refer to elements of the point cloud and never to vertices. Next we give 
definitions and some properties of barycentric coordinates and Bernstein polynomials. 


More information and basic proofs can be found in Farin [11]. 


Definition 1.1. Given a simplex A = (vo, v1,..., Va) the barycentric coordinates of a 


d 

point p € R@ are (Ao, A1,..., Aa) such that p = Agvo +Arvi +++ -+Aagva and S> A; = 1. 
i=0 

Although an abuse of notation, we write A,(A, p) as the barycentric coordinate of p 


in the simplex A with respect to the vertex v. 


The usefulness of the barycentric coordinates is based on their simplicity and 
flexibility in any d. Each barycentric coordinate of p is essentially the weight of the 
corresponding vertex for p. Furthermore, the barycentric coordinates with respect to 
a nondegenerate simplex are unique for each point because we enforce that > be | 
in the definition. When p is near a vertex we have that the corresponding taco 
coordinate will be large, and p far away from a vertex will have a small barycentric 
coordinate. In fact, the barycentric coordinates of the vertices of a non-degenerate 
sitriplex are Ave(A,.Vo) = "(Oj 850) Ag (Agvn) =: (05.1, Oy cece oy Ag (Ay va) = 
(0,0,...,1). If all of the barycentric coordinates are positive, then p lies inside the 
convex hull of the vertices. On the other hand, p lies outside when at least one 
coordinate is negative. Furthermore, barycentric coordinates are affine invariant. 
That is to say that they do not change when the simplex and point undergo the same 


affine transformation. 


The barycentric coordinates A; to Aq are found by solving the following system: 


ry 
A2 
(v1 -vo V2—-—Vo -c° vi=ve) : = Pp. 
Xa 
The matrix is invertible since the v; are the vertices of a non-degenerate simplex. 
d 
Then the Ag is found by Ay = 1 — y, A;. For a triangle (vo, V1, V2), the formulas are 
i=1 
simply 
area(p, V1, V area(Vo, Pp, V area(Vo, V1, 
i ee Se DD im 
area(Vo, Vi, V2) area(Vo, V1, V2) area(Vo, V1, V2) 


Barycentric coordinates are used to define the Bernstein polynomials, which are 
utilized in Chapter 3. The most important property of the degree n Bernstein poly- 


nomials is that they form a basis for all the polynomials of degree n. 


Definition 1.2. Let (Xo, \1,..., Aa) be the barycentric coordinates of p. The Bern- 


stein polynomials of degree n over a simplex of degree d are defined by 
Bi(p) = (") APAE AY 
g n n) 

where t = (to,..,ta) € N4*?, Sot; =n, and : 


Example 1.3. The Bernstein polynomials of degree n over the interval {a, b] are: 


Fe (eae 


t 


Some properties of the barycentric coordinates directly carry over to the Bernstein 


polynomials: 
e affine invariance 
© Bip) = Oot +20)" =1 
e p in the simplex if and only if Bf (p) is positive for all t. 
All of the Bernstein polynomials are symmetric with respect to the center of the 


simplex. 


1.3. TRIANGULATIONS AND NEWEST VERTEX BISECTION 


This section explains the method we use to partition the domain X containing the 


point cloud D in R?. The first topic is to generate partitions using the newest vertex 
bisection method of triangle subdivision. We discuss the natural master tree that 
corresponds to the partitions generated by the newest vertex bisection. ‘Together, 


the tree structure and the subdivision method allow for multiresolution analysis of 
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Figure 1.1 The domain X with initial partition Po. 


the subdivision surfaces. Then the topic shifts to adaptively choosing a particular 
partition from the newest vertex bisection method based upon the aggregate quantities 
of the point cloud. The last topic discusses the usefulness of conformal partitions and 


the procedure to ensure that the partition is conformal. 


Definition 1.4. A partition P of X C R@ is conformal if any two simplices in P are 


either disjoint or intersect at a simplex of lower dimension. 


Newest Vertex Bisection 


Newest vertex bisection has been used successfully in adaptive finite element meth- 


ods (AFEM) to establish convergence and rates of convergence (see, Morin [13] and 


[5]). Furthermore, newest vertex bisection has generalizations to R® by tetrahedral 


bisection [1] and to R@ by simplex bisection [14]. The tree structure benefits carry to 


R? in [6]. 


For our convenience we assume X is a square, or in R? a d-dimensional cube. The 


newest vertex bisection method subdivides a triangle by bisecting the edge opposite 
the newest vertex. So to start, we need an initial partition Pp and an assignment of 
newest vertices. To keep the partition simple, we define the initial partition of the 


square X to be two right-isosceles triangles as seen in the figure. 
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VB 


VL VR 


Figure 1.2 The reference triangle with vertices vp = (0,1), vz = (—1,0) and 
VR>= Ch 0). 


VL Vy VR 


Figure 1.3 Subdivision of the triangle A = (vg, vz, vp) into its children 
Ao = (Vs, Vp, Vz) and Ay = (vx, VR, VB) by newest verter bisection. 


We label each vertex corresponding to a right angle as the newest vertex in its 
respective triangle. Newest vertex bisection with these initial conditions gives us that 
all triangles in the partitions will be right-isosceles. Furthermore, the newest vertex 
of a triangle will always be at the right angle. 

Let A be a right-isosceles triangle. Then A is uniquely identified by its vertices 
(VB,VL,VR), where vg is the newest vertex. Any right-isosceles triangle can be 
transformed to the reference triangle with vertices at (0,1), (—1,0), and (1,0). The 
vertex corresponding to (0, 1) is labelled vg , (—1, 0) is labeled vz, and (1, 0) is labeled 
VR. 

We subdivide the parent triangle A = (vg, vz, vr) by bisecting the edge opposite 
the newest vertex and adding a new vertex v, at the midpoint. The left child is 
Ao = (Vs,VgB,Vz) and the right child is A; = (v,,vR,vg). The vertex v, is the 


newest vertex of the children, as the name would suggest. Bisecting the hypotenuse 
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of a right-isosceles triangle always results in a right angle and the finer triangles 
created are right-isosceles. Therefore the newest vertex is always at the right angle 
of a triangle and all of the triangles resulting from further subdivision will also be 
right-isosceles. If a different initial partition or labeling was chosen, then we may not 
have that all of the triangles are right-isosceles. The framework shown here will be 
relevant for more general partitions when the initial setup is different. 

Now we will describe the tree structure that accompanies the newest vertex bisec- 
tion. The initial partition P) consists of two triangles Ap and A,. In fact, we write 
Po = {Ao, Ai}. Construction of our master tree T begins with an empty root Ag. 
Then the left child of the root is Ap and the right child is A, which are the triangles 
that make up Pp. Next, Ag is subdivided into its children Agg and Ap; by the method 
described above. The master tree is expanded by appending Agg and Ag; to Ag in T. 
Continuing to subdivide subsequent generations by newest vertex bisection appends 
more and more children nodes to 7’. Thus, Tis an infinite rooted tree that is strictly 
binary. The set of nodes of T are all of the triangles which can be obtained by a 
sequence of subdivisions starting from Po. Once we have Po and the initial newest 
vertex assignment, then the master tree 7’ is predetermined. 

Full binary trees are defined to be binary trees in which all nodes have either two 
children or no children. We define a pruned subtree to to be a subtree that contains 
the root. Let S be the set of all finite, full, and pruned subtrees of 7’, and let P be 
the set of all partitions of X generated by the newest vertex bisection. For any tree 
S let N(S') be the set of the nodes of S and £(S) denote the leaves of S. 
Definition 1.5. L C N(T) is a partition in P if and only if |) A; =X and the 


A;EL 
intersection of any two distinct triangles in L has planar Lebesgue measure zero. 


Lemma 1.6. Define R: S > P by R(S) = P whenever L(S) = P. Then R is 


bijective. 


13 


Proof. Let P be any partition in P. Then by definition P C N(T). Create S, which 
will be the pre-image of P, by appending to P all of the ancestors of the triangles 
in P, so that S is a pruned subtree of 7’. If any triangle in P were the parent of 
another triangle in P, then P would not be a partition. So P = £(S). It remains to 
show that S is full. Suppose S were not full. Then there exists a node A € S with 
children Ap and A, in T such that only one child is in S. Without loss of generality, 
assume Ao is in S. Then A, Aj, and the children of A; are not in P. This implies 
that P does not cover all of X. Therefore S € S which implies R(S) = P. Hence R 
is surjective. Let S, and Sj be in S such that R(S,) = R(S2). Then £(S,) = LS). 


Since they are full binary trees with the same root and same leaves, we have that 


S, = So. So R is bijective. 


Now that the one-to-one correspondence between S and P has been established, 
we introduce two operations that will allow us to describe families of partitions: 
splitting and merging. Let A be a triangle with children A; and Ay. For S € S 
and Aj, A, € L(S), the operation merge (Aj, Ay) removes both A; and Ag from 
S, which implies that A is now in £(9). For S € S and A € L(S), the operation 
split (or subdivide) A, does the opposite in that it appends the children of A in 
the master tree to S. Merging two triangles of a partition decreases the number of 
triangles in the partition by one, while splitting a triangle increases the number of 
triangles by one. In fact, any partition in P can be obtained from any other partition 
in P by a sequence of splits and merges. 

The structure of the master tree and the subdivision method allow for multires- 
olution analysis, since the partitions of finer levels maintain the vertex locations of 
the coarser levels. Subdivision by newest vertex bisection starts with a base mesh of 
two triangles and recursively obtains finer meshes, not necessarily uniformly. Finer 
meshes correlate to increased depth in the master tree. Each time a triangle is sub- 


divided, the location of the vertices in the parent remains fixed and a new vertex is 


14 


added, which implies that a fine mesh contains all of the coarser meshes. A surface can 
be described at varying depths of the master tree, specifically coarse approximation 
to fine. Applications of multiresolution meshes include progressive transmission and 
level of detail control. Progressive transmission begins with a low resolution approxi- 
mation with very few details, and then progressively improves by adding more details. 
Level of detail is controlled by the addition or removal of subdivision branches. 

Of primary interest to any implementation of this newest vertex bisection method 
is having an efficient way to calculate and maintain the quantities we store in each 
node of the tree. The simplest way is also the most computationally expensive: when 
a triangle is split, sort the points from the parent into each child and then calculate 
the quantities for each child directly from those points. A more efficient algorithm 
in which a tree skeleton is constructed and the quantities are calculated directly on 
the finest level and then are added together for coarser levels is described in [6]. 
Additivity of the quantities is a necessary property in order to use this more efficient 
algorithm. The algorithm in Chapter 3 needs neither a tree skeleton nor continual 
sorting of the points. We calculate all the quantities of the finer level directly from 


the coarser level. 


Adaptive Partitions 


Using the information from the point cloud, we adaptively choose a partition P4 € P 
that satisfies the aim of the application. That is: we tailor the partition to both the 
specific requirements of the application and the information present in the point cloud. 
For instance, the partition should be subdivided more (making a finer mesh) near ar- 
eas of rapid change or important details and subdivided less (remaining coarser) over 
areas of slight change. The criteria driving the choice of P, are motivated by the fac- 
tors that will produce a quality approximation of s,. There are two main approaches 


to arrive at an adaptive partition: Coarse-to-Fine and Fine-to-Coarse. The Coarse- 
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to-Fine approach is driven by the subdivision criteria discussed in Section 2.3 while 
the Fine-to-Coarse approach is driven by the merging criteria. The tree structure al- 
lows for a combination of the two approaches using both splitting and merging. The 
adaptive partition P, is the partition resulting from repeatedly subdividing Po until 
none of the leaf nodes meet the subdivision criteria, or similarly for Fine-to-Coarse: 
repeatedly merging a finest-level-of-detail partition until none of the leaf nodes meet 
the merging criteria. 

In the Coarse-to-Fine approach, we start with the initial partition Pp and, based on 
subdivision criteria derived from the quantities of the triangles, we mark triangles for 
subdivision. Then all marked triangles are subdivided. When conformal partitions are 
desired, we also split neighboring triangles as needed to resolve hanging vertices, even 
when those triangles do not meet the subdivision criteria themselves. See Conformal 
Partitions in the following section. After each stage of subdivision, the quantities 
are recalculated. Triangles are marked and subdivided repeatedly until all of the 
triangles in the partition satisfy the subdivision criteria. One benefit of the Coarse- 
to-Fine approach is that only the necessary triangles in T’ need to be calculated. For 
very large point clouds, this saves a significant amount of computation. 

In the Fine-to-Coarse approach, we start with some large finite subtree of the 
newest vertex bisection tree S C T, in which S is large enough that no finer details 
are necessary. Then the quantities are calculated for the leaves of S. Triangles 
are then merged based on the merging criteria. If conformal partitions are desired, 
then only diamonds and boundary triangles are merged, and forced mergings may 
occur. The Fine-to-Coarse approach allows for easy encoding and computation of the 
quantities. 

The criteria used to mark a triangle for subdivision or merging vary and can be 
tailored to the specific goals of the application. Most commonly an error term is 


defined for a triangle. The error could be the variance, the sum of the residuals, 
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and so on. Whatever criteria is used, it is important that the quantities involved are 


easily computed. 


Conformal Partitions 


Applications in which the sensed surface is supposed to be continuous desire that 
the surface approximation is also continuous. The partitions generated so far might 
not produce a continuous approximation. More specifically if a vertex of one triangle 
intersects the edge of another triangle, then a gap may occur between the vertex 
and the edge. Such a vertex v is called a hanging vertex of P, as v appears in 
the interior of an edge in P. In conformal partitions (Def. 1.4) hanging vertices 


are not permitted, which guarantees that the approximations built on conformal 


partitions are continuous. In R?, triangles of a conformal partition P are permitted 


to intersect only at a vertex or an edge. This section discusses two procedures of 
receiving conformal partitions: completion after subdivision and force-splitting. If 
the partitions are not required to be conformal then it is not necessary to implement 
either procedure. In Chapter 2, we apply the force-splitting procedure to produce 
a conformal partition. However, our sensed surface is composed of different objects 
(e.g., tree, building) meaning only parts of the surface would be continuous. So we 
adjust for discontinuities using criteria discussed in Section 2.1, and our final partition 
will not be conformal. 

Completion after subdivision is a procedure to turn a non-conformal partition into 
a conformal one. This is a typical procedure used in adaptive finite element methods 
that mark triangles for subdivision and then split all marked triangles. After all of 
the subdivisions are made to satisfy the subdivision criteria, the method finishes with 
a completion step. In the completion step, the partition is further subdivided until 
there are no hanging vertices. The number of triangles subdivided in the completion 


step is shown to be bounded by a constant multiple of the number of triangles in the 
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uncompleted partition (see, [5]). 

The second method consistently maintains the conformality of the partitions using 
force-splits (see, {10]). Thus force-splitting is particularly useful for the Coarse-to- 
Fine approach to adaptive partitions. In order to prevent hanging vertices along the 
way, only certain triangles are eligible for immediate subdivision, and when a triangle 
is not eligible a neighbor will be forced to subdivide as well. A benefit is that the 
partition is conformal after each subdivision. For the convenience of the reader, we 
include some details here. 

We remain in the same setting as in Section 1.3 using newest vertex bisection. 
Recall that the cells are right-isosceles triangles and are marked for subdivision based 
on the subdivision criteria in Section 2.3. Now we force additional triangles to split. 
Although we desire to subdivide a triangle A, we may have to subdivide other trian- 
gles first so that hanging vertices are avoided. A triangle A is eligible for subdivision 
under one of the two following circumstances. First, A is a boundary triangle, which 
is defined as a triangle whose hypotenuse is on the boundary of X. Second, A shares 
a hypotenuse with another triangle from the same level in 7. This case is referred 
to as a diamond, consisting of A and the triangle sharing the same hypotenuse. In 
the case of a diamond, the subdivision of A forces the other triangle to immediately 
split. 

The idea is that in order to subdivide a non-eligible triangle A,, we first must 
subdivide other eligible triangles until A, becomes eligible. Let A, be a non-eligible 
triangle marked for subdivision. Then we define the following recursive procedure 
terminating when A, becomes eligible and the diamond is split. Since A, is not 
eligible we have that A, has a neighbor A’, from a coarser level. Geometrically, the 
depth of A, is required to be exactly one more than the depth of A‘. If AY is eligible 
for subdivision, then we subdivide A’, making A, eligible for subdivision. It is possible 


that A‘ is not eligible for subdivision. In this case we mark A’ for subdivision, and 
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recursively repeat the process above marking until we arrive at a triangle that is 
both marked and eligible. The triangles encountered during this recursion continue 


to decrease in depth, so the recursion always terminates. 


Singular Value Decomposition 


A singular value decomposition (SVD) provides a convenient way for breaking a 
matrix into a product of simpler matrices. The SVD factors any m x n matrix A into 
the form UNV? where U and V are unitary and D is diagonal. The diagonal entries 


of © are called singular values and are often denoted by o; for i = 1,..,n. 


R2x2 


To visually understand SVD, consider the geometry in R? for A € a full rank 


matrix. The image of the unit circle under multiplication by A is an ellipse. Then 
the singular values of A are the lengths of the semi-principal axes of the ellipse. The 
matrices U and V can be thought of as rotation matrices. It is possible for some 
singular values to be zero or near-zero. In our geometric perspective, a near-zero 
singular value means that the minor axis of the ellipse is so small that the ellipse is 
flattened very close to a line. Thus, in general dimensions, having a zero singular 
value implies that the data in A lies in a space with dimension less than n. 

For our purposes, the matrices with singular values near or equal to zero are of 


particular interest. In these circumstances one can consider the truncated SVD. Let 


t be a given threshold for the singular values, and A € R™*”". Then suppose A has k 


singular values greater than t. Then only the k column vectors of U and the k row 
vectors of V? corresponding to the k largest singular values of © are calculated. The 
rest of the matrix is discarded. Thus A = A = U;,5,V,’. The approximate matrix 
A is the closest approximation to A that can be achieved by a matrix of rank k. 
Naturally, if k = rank(A), then A = A. 

One application of SVD is solving a system of equations. Assume Ax = b and we 


are solving for x. Then the SVD of A gives us that there exists a diagonal matrix ¥ 
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and two unitary matrices U and V such that A = UNV". (Side note: When A is real 


symmetric, we have that U = V = U7 and that the entries of © are nonnegative.) 


Ax =b 
UXV'x=b 
(V"x) =U*b. 


Assuming the singular values of A are nonzero, solving for V’x is simple because © is 
diagonal, and let S denote this solution. Using the truncated SVD would also avoid 


any issues with solving for Vx. Then 


Vix=S > x=VS. 


There are several algorithms for calculating the SVD with stable performance for 
moderate dimensional matrices. We will be using the SVD only in three dimensions, 
so there is no problem using the SVD standard in MATLAB or other computation 


software. 


Hausdorff Distance 


The Hausdorff distance measures how far two subsets of a metric space are from each 
other, and is most often used when comparing two representations of the same three 
dimensional object. We use the Hausdorff distance when comparing our point cloud 


approximation to the original surface. 


Definition 1.7. Let X and Y be two non-empty subsets of a metric space (M,d). We 


define their Hausdorff distance dy(X, Y) = max {sue inf d(x, y), sup inf d(z, i}. 
rEX yeY yeY rEX 


For our purposes M will be R® and d is the Euclidean distance. Simply put, the 


Hausdorff distance is the maximum of the distance of the farthest point in X to Y 
and the distance of the farthest point in Y to X. Note that the Hausdorff distance 


is not generated by a norm. 
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sup inf d(x, y) 
reX ¥EY 


d(x, y) 


sup inf 
yeY rex 


Figure 1.4 Illustration of the Hausdorff distance between the sets X and Y. 


Often the L?,L', or L® errors are used, but we argue that the Hausdorff distance 
offers a better measurement of the similarities between two surfaces in terrain approx- 
imation modeling. Take for instance two scenes that differ only by a telephone pole. 
A telephone pole is used in the illustration because the base of the pole is very small, 
and the height is large. The L? and L' errors would be small, yet the Hausdorff would 
be large. The L° error would also be large; however, a slight horizontal shift in the 
approximation could result in very large L° error. Therefore, we use the Hausdorff 
because differences like telephone poles are significant in our applications. 

The presence of outliers greatly impacts the Hausdorff distance. Outliers may be 
present in the point clouds we analyze, but we use the Hausdorff metric to compare 
our surface approximation to the actual surface. In this case our approximation has 
already decreased the influence of the outliers, so the comparison of the surfaces are 


fairly resistant outliers. 


1.4. SUBDIVISION 


Subdivision surfaces are becoming the standard surface models, and are used in a wide 
range of applications. Subdivision is easy to implement and computationally efficient. 
In this section we discuss splines [8] as a one dimensional example of subdivision and 


4-8 subdivision [15] as a two dimensional example. Then we explore the idea of 
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polynomial reproduction and the influence it has on the development of subdivision 


schemes. 


Splines 


Splines are smooth piecewise polynomial curves of some designated degree with high 
continuity at the junctions. We show here how to define splines and how they can 
be generated through subdivision. We start at the very beginning with B-splines. 


B-splines are defined through repeated convolution, 


t Oot 
Ba) = (1.4) 


0 otherwise 


and then 
B(x) = BY(x) @ B(x) = i BY(t)B™ “(x — t)dt. (1.5) 


From this definition there are some important properties that follow. The first is that 
the B-spline of degree n is C"~! continuous, and is a direct consequence of convolution 
with B!. The second property is that B-splines obey a refinement equation. That 
is, a B-spline can be written as a linear combination of dilations and translations of 


itself. The refinement equation is given by 


= ye (7) Br ee— (1.6) 


+e (ren (res 


For instance 


(B? (2x) + 2B? (2x — 1) + B?(2x — 2)). 


The spline curve y of degree m is written as a linear combination of shifted B-splines, 


and we call each coefficient a control point. Then 


Ze 


Notice that each control point p; only influences t € [i,i-++m]. To simplify our notation 
let p be the column vector of control points and let B(t) be the row vector of the 
translates of B”. Thus y is written y(t) = B(t)p. The refinement equation (1.6) 


becomes B(t) = B(2t)S, where S is a matrix with entries 


1 m 
Saitkk = sal"). 


This is refinement of the basis elements. Therefore, we can write y as a linear com- 


bination of B-splines B(2t) and new control points Sp 
y(t) = B(t)p = B(2t)Sp. 


This shows that basis refinement corresponds to control point refinement. Repeating 


the refinement gives us 


y(t) = B(t)p® (1.7) 
= B(2t)p' = B(2t)Sp° (1.8) 
(1.9) 

= B(2/t)p’ = B(2’t)S’p° (1.10) 


Thus the relationship given between control points at level j of subdivision is p/+! = 
Sp’. Now consider the control points themselves. The graph of the control points with 
lines connecting consecutive points is a piecewise linear approximation of the curve 
y. The initial control points p are a coarse approximation, and each multiplication 
by S gives us a finer approximation. For the purposes of computer visualization, one 
draws the piecewise linear curve to the desired level of detail, instead of drawing the 


curve ¥. 


4-8 Subdivision 


Subdivision schemes on surfaces are defined depending on the rule of partitioning and 


the goal is to ensure higher smoothness of the resulting surface. Usually, this is done 
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through considering well-known smooth bell-shaped functions y to be a generating 
function for the regular schemes, meaning that the limiting surface y will have the 
representation 
yz) = 2 Paolo —k), 

where px are the initial values at the control points. Often the generating function of 
choice is a box spline. For instance, three-directional quartic box splines lead to the 
Loop subdivision scheme [12], and the tensor product biquadratic and bicubic splines 
lead to Doo-Sabin [9] and Catmull-Clark [7] subdivision respectively. The box splines 


are defined similar to the B-splines above. Let [d)...d,,] be a set of directions. The 


box spline B™(x),x € R? can be computed using the following recurrence: 


Bi(x) = [ Brix -tajat, 


with B° being the delta function. The directions can be repeated and usually three 
or four directions are chosen. The directions for four-directional box splines are (1, 0), 
(0,1), (1,1), and (1,-—1). The 4-8 subdivision scheme uses bisection refinement as 
the one in newest vertex bisection defined in Section 1.3 but with the additional 
requirement: if an internal edge is bisected, then both its adjacent triangles are 
bisected. For the unit square the usual initial partition consists of four triangles 
received by adding both diagonal of the square and declaring their intersection point 
to be the newest vertex for all four triangles (it can be received from the initial 
partition in Figure 1.1 by subdividing its internal edge). On each step we subdivide 
all the triangles to receive the next partition. The signature property of this partition 
is that the valence (number of edges meeting at a vertex) of the internal vertices is 
either 4 in case of newest vertices, or 8 otherwise. This is the reason for the name “4-8 
subdivision scheme”. Full description of the 4-8 subdivision scheme is given in [15]. 
Here we show only the schemes for the regular internal vertices noting that there are 


special rules for treating the boundary vertices, as well as the cases of extraordinary 
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Figure 1.5 Mask for the addition of vertex to a face in the 4-8 Subdivision Scheme. 
l l 


N71. ZN) 
TAN) N77 


Figure 1.6 Masks for the regular case 4-8 Subdivision Scheme when updating a 
vertex. 


vertices which can be introduced in the initial partition with valences different from 
4 or 8. The mask for calculating the value at the new vertex at the middle of the 
diagonal is given in Figure 1.5. The masks for updating the values at the existing 
verices are given in Figure 1.6. Note that in both figures the partitions are shown 
before the new bisections are made. 

The generating function for the 4-8 subdivision scheme on a regular partition is 
the four-directional box spline for which each of the four directions is taken twice. It 
is a piecewise polynomial function of degree 6 and has continuous derivatives of order 
4. Furthermore, the subdivision surfaces produced by the scheme are C* continuous 


almost everywhere, except at extraordinary vertices where they are C' continuous. 


25 


CHAPTER 2 


ASSIMILATION AND PROCESSING OF POINT CLOUDS 


This chapter builds on methods from Learning Theory to develop techniques and 
algorithms for handling the specific and difficult case of irregular point clouds. The 
main idea is to take methods that have have proven to be successful for relatively 
regular point clouds and extend them to handle irregular point clouds. We preserve 
the efficacy of these methods by transforming our point cloud to a different topology, 
one in which the point cloud is semi regular. In particular, the point clouds received 
from stationary LiDAR are more regular when represented in spherical coordinates. 
Even in spherical coordinates, there still could be some gaps in the data where the 
LiDAR did not receive a return. The sensor location is the origin, and the points 
are given in angle-angle-distance. Indeed, most LiDAR hardware naturally outputs 
the readings in spherical coordinates, only to later be transformed to Cartesian co- 
ordinates to assimilate with existing data sets. Here we only consider point clouds 
resulting from a LiDAR scans, and do not consider dense and regular (in Cartesian 
space) point clouds resulting from several scans or perfect knowledge. One area of 
our investigation is to assess how the Cartesian and spherical coordinate systems in- 
teract. Another area is how to diminish the effects of missing data and jumps in the 


distances from the sensor. 


2.1 THEORETICAL FRAMEWORK 


The practical algorithms developed in this chapter realize ideas expressed in the 


work of Binev, Cohen, Dahmen, DeVore, and Temlyakov in Universal Algorithms 
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for Learning Theory. Here we discuss the main results for approximations using 
piecewise constants [4] and piecewise polynomials [2] over adaptive partitions. Let 
D = {(x,y)}™, be the point cloud. The convergence of the approximation 8 
to the regression function s, is usually evaluated in probability or in expectation. 
The approximation § converges to s, in probability if for all 7 > 0 the limit of 


Prob{||s, — 5|| > 7} is zero as the sample size approaches infinity. When the limit 


of E( 


S$» — 4||?) is zero as the sample size approaches infinity, the approximation is 


said to converge in expectation. Given a partition P of X, the piecewise constant 
approximation § is the solution to the least-squares problem. More specifically, § is 


the piecewise constant function over P that is the minimizer of 


which implies that,over each cell of P, § is simply the empirical average. ‘Their 


algorithm consists of the following steps: 
e Compute the error over each cell. 


e Threshold at level tr, := &,/ wen to obtain the smallest proper subtree of the 


master tree containing cells with error larger than t,, (& a constant). 
e Complete the subtree and obtain the corresponding partition P. 
e Compute s by empirical risk minimization over P as described above. 


Technically, the regression function s, needs to be in A7()B’. The approximation 
class A? can be viewed as a smoothness space of order y > 0 with smoothness mea- 
sured with respect to px. Given v > 0, B” is a smoothness space which measures the 
regularity of a function f by the size of the adaptive partition. Details for the algo- 
rithm and the smoothness spaces can be found in [4]. The main result for piecewise 


constant approximation is the following theorem. 


yas 


Theorem 2.1. Piecewise Constant Approximation [4] Let 3,7 > 0 be arbitrary. 
Then, there exists ko = Ko(B,y,M) such that if & > ko, then whenever s, € AT) B” 


for some v > 0, the following concentration estimate holds 


l BI 
Prob {ts —al|>é ( =”) < Cm? (2.1) 
as well as the following expectation bound 
2 sl?) < c ( 8™ oe 2.2 
(IIs. — 3?) < C | —— (2.2) 


where ¢ and C are independent of m. 


The order of the approximation is optimal save for the logarithmic term. The 
analog of Theorem 2.1, when replacing piecewise constant with piecewise polynomial, 
holds with the addition of some regularity conditions on px. The counterexample, 
given in Section 3 of [2], shows that the analog is not true without the regularity con- 
ditions. The regularity conditions imposed upon the measure are easily satisfied for 
measures equivalent to Lebesgue measure, but not in the case of gaps in the measure. 
Details of the conditions can be found in [2]. When px is a very irregular measure, 
the theorem does not imply optimal results for piecewise polynomial approximation. 
Regardless of the regularity, the algorithm still yields an approximation §. 

The empirical risk minimization does not perform well in probability for the coun- 
terexample in [2]. Consider px to be an irregular measure with high concentration on 
one side of the subdomain, and low concentration on the other. Given one realization 
of the point cloud we may only have points on the side with high concentration of mea- 
sure and no points on the other. However, in another realization the point cloud may 
have a few points on the side with low concentration. In both instances the algorithm 
would approximate the surface, but would receive significantly different results. This 
is the essence of the counterexample showing the instability of the approximation. 


Our practical algorithm uses two techniques to suppress the instability. 
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1. Truncated SVD: The truncated singular value decomposition (SVD) is used in 
the vertex value calculation because it automatically approximates over both 
uniformly distributed data and data lying on a hyperplane of lower dimension 
(see, Section 1.3). Essentially, this finds a subspace on which the regularity 


conditions hold, and ignores the part that might cause problems. 


2. Special Subdivision Criterion: We impose a subdivision criterion that splits a 
triangle if the distance between the center of mass of the subdomain and the 
center of mass of the points in the subdomain is relatively large (see, Section 
2.3). This criterion reduces the number of triangles in the adaptive partition 


that may have this issue. 


Notice that the first technique modifies the vertex approximation, and the second 
technique modifies the adaptive partition. 

Often the surface represented by the LiDAR point cloud is discontinuous. The 
discontinuities arise whenever one object ends, from the vantage point of the sensor. 
Examples are the corner of a building or the outline of a tree or telephone pole. Over 
these discontinuities, large jumps in the distance from the sensor occur over small 
changes in the angle. If modeled by a continuous surface, then there would be long 
skinny triangles along the line-of-sight of the sensor. This hinders the quality of our 
approximation, so we relax the conformity and continuity requirements in two ways. 
The first is that we remove empty triangles from the partition. Where there are no 
points, there is no information from which to construct the surface. The second is 


that we permit a vertex to have more than one value near discontinuities. 


2.2 SPHERICAL LEARNING SETTING AND NOTATION 


This section sets the stage of both the spherical and Cartesian coordinate systems in 


which we process point clouds. The LiDAR sensor location is fixed and rotates to scan 
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over the range —7/2 to 7/2 latitudinally and —7z to m longitudinally. Note that some 


devices are not capable of rotating the full range. Let R? = X xY refer to the standard 
Cartesian coordinate system with 71,72 and y. Define S = A x R to be the spherical 


coordinate system with 0,¢, and r where A = {(6,¢) : 6 € [—7, a], ¢ € [—7/2, 7/2]} 


and R Cc R. We have that —7 < 6 < 7 is the longitude, —7/2 < @ < 7/2 is the 


latitude, and r > 0 is the distance from the origin, i.e. sensor location. Recall the 


Cartesian coordinates for a point p are of the form (x,y) with x € R? andy € R. 


Recall from Preliminaries 1.1 x is referred to as the location and p = (x, y) is a point. 
In spherical coordinates p = (a,r) where a = (6,¢) is the direction from the origin 
(sensor) and r = r(a) is the distance to the point from the origin in the direction a 
and is considered to be a function of # and ¢. The angle @ is the angle of rotation 
about the vertical axis (y direction) from the positive x, axis. The angle ¢ is the 


angle of elevation from the x; and x2 plane. The standard conversion formulas from 


S to R® are: 


x1 = 7 cos(¢) cos(@) 
xr2 = rcos(¢) sin(@) (2.3) 


y =rsin(¢) 


and from R? to S: 


= eae 
P y 
@ = arcsin (2) (2.4) 


6 = atan2(z2, 71) 
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More specifically 
arctan (12/21) x1 >0 
arctan(%e/11) +7 21 <0,%2>0 


arctan(%q/r1) —m 21 >0,2%2 <0 
atan2(x2, 21) = 


m/2 = O50) 
—1/2 rv, = 0,22 < 0 
undefined t7, = Ost = 0 


is a common function used in computational science to find the angle using inverse 
tangent and taking into account the quadrant the angle lies in. 

We assume that D consists of independent random observations that are identi- 
cally distributed according to p on A x R, and the points in D can be represented in 
either coordinate system. In this case the marginal probability measure is denoted 


pa and the regression function is 


sp(a) = | rdp(rla). 


Here we are approximating the function s, = s,(a). For a spherical triangle Ag C A 


the quantities that drive our approximation are the following: 
e Spherical Moments: 
| ldpa, | dpa, | odpa, | rdpa, and | Ir|/?dpa 
As As As As As 
e Cartesian Moments: 
| nidp, | radp, | neadp, | vidp, | a3dp 
AsxR AsxR AsxR AsxR AsxR 


d sh d if do and i: 2d 
cet I acne? p ed ey hile p 


e Extrema: 


max{r:(a,r)€ Ag} and min{r: (a,r) € Ag} 


dl 


The spherical and Cartesian moments have the desired property of additivity. The 
extrema are not additive, yet are easily determined. When merging two triangles, 
the maximum r-value of the parent is the maximum of the two children maxima, and 


the same for the minimum. 


2.3. SURFACE APPROXIMATION 


Adaptive Partition 


The first step in building our surface approximation is to partition the domain. To 
arrive at our adaptive partition P, we follow the procedures outlined in the prelim- 
inaries. Since the point cloud is regular in spherical coordinates, we partition the 
S domain, A, which is a square encompassing the point cloud, using newest vertex 
bisection with force-splitting as explained in Section 1.3 to some finest level. The 
triangles forming a partition in A correspond to spherical triangles in R® where the 
edges are partial great circles. 

There is one difference to note between partitioning A, and the Cartesian domain 
X. In the preliminaries, the domain X has clearly defined boundaries, and outside of 
those boundaries there are no connections. However the triangles on the boundaries 
of A with |0| = 7 are neighbors and share vertices or edges. This requires the 
subdivision scheme to link these triangles and include the connections for the purpose 
of force-splitting. The connections are also needed for any calculation based on the 
neighborhood around a vertex. 

The second step sorts the points of D into the triangles of the finest level based 
on their 6 and @¢ coordinates. Then the aggregate quantities are calculated on the 
finest level directly from the points. After this, the algorithm never uses the points 
themselves. Working back up the master tree, the additive quantities of the children 


are added together to give the corresponding quantities of the parent. For the max- 
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imum and minimum, the extrema of the parent is the most extreme extrema of the 
children. Next we give the criteria which adaptively determines a partition P, € P. 

The main motivation in deciding our criteria is that we want to have a stable 
and accurate approximation using the easily computed quantities. Our subdivision 
criteria marks triangles for subdivision if the variance is large or the center of mass 
of the triangle is far away from the center of mass of the points it contains. Let Ag = 
(Vo, Vi, V2) be a spherical triangle whose vertices have S coordinates (0%', é%',r%). 
Furthermore let {p}”_, be the points in Ag, with S coordinates (0, 6,r). The 
average r-values of the points in Ag is denoted by 7. Similarly, @ and ¢ denote the 


averages of 6 and ¢. Then the variance of the r-values of the points in Ag is 


= Jas Ir Pdpa _ G a ey 
Jas ldpa Jas ldpa 


For the purposes of the subdivision criteria, the center of mass of the points in As 
Qvo + gv + v2 gre + ris + pr? 
, . Define 


Var,(As) 


is (@,@) and the center of mass of Ag is ( 5 3 


d(As) as the Euclidean distance in S between 


gvo 4 v1 +4 Qv2 pre + o™} 1) 
3 ; 3 ; 


(0,6) and ( 
Therefore, a triangle As is marked for subdivision if either Var,(Ag) or d(As) is larger 


than some respective threshold. We list the specific criteria and thresholds used in 


the numerical results in Section 2.5. 


Vertex Approximation 


In this section we calculate the r-values for the vertices in the adaptive partition P, 
generated in Section 2.3 using the aggregate quantities of the points in the triangles. 
First we calculate all the local r-values for a vertex, where each local value uses only 
the quantities from one triangle containing the vertex. Then we combine the local 
r-values into a single global r-value. Under certain circumstances, we may assign two 


distinct r-values to a single vertex. 
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Local Approximation 


The idea of this section is to give an algorithm to approximate the value of a vertex 
v over a single subdomain As that contains v. We use the method of least squares 
to fit a plane over the triangle Ag based on the aggregate quantities. Higher degree 
polynomials are not used because their calculation requires denser point clouds than 


what we receive from real world LiDAR data subjected to a relatively low sampling 


rate. Therefore we transition from S to the real space R®. We cannot simply express 


the vertices in R?, as they are in the form (0, ¢) with unknown r. Until the r-values 
are found, each vertex v could be thought of as a ray emanating from the origin 


and denoted by v. First we find the plane of best fit, A, through the points in 


As = (Vo, Vi, V2) using R® coordinates. Then the intersection of K and v; provides 


r; for each i = 0,1,2. In the case that K and v; are almost parallel, the r-values lie 
outside of the range of the r-vales of the points in Ag. If so, we clamp the r-values 
to the maximum or minimum. 

Let F = {(x,y) € D: x € As} and n denote the cardinality of F. Then we have 
the local points F = {x’, y'}",. Specifically, we let E = {x’}"_, be the location of the 
points in F’. We can treat the data (y) as both a vector of values, y = [y’,..., y”], and 
as a real valued function of the discrete domain FE, where y(x’) = y’ for i = 1,...,n 


Let f and g be two real valued discrete functions defined on E’. Then we define the 


discrete inner product in the usual way 


Fe = Ta Fle > Fee (2.5) 


xe i=l 


This discrete inner product relates to the following inner product for continuous 


functions f and g over Ag in S. 


1 
(f, 9a, = Tit bape £900 


The inner products are calculated from the Cartesian moments for Ag. We define 


our linear least squares method that approximates s,(v) using inner products. 
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Here we describe the least squares setup. Let K(F’,x) be the linear approximation 
of at x. Traditionally, K is of the form 
d 
KEXP, Rep Ci, 28,64) = 2, C55(%) 
= 
where the c; are d+ 1 adjustable parameters of AK, and the @; are a basis for linear 
functions over the domain. Start by letting @9(x) = 1. Then the remaining linear 
functions ¢; (j = 1,...,d) are chosen to be orthogonal to ¢. So we set ¢,(x*) = 
(xf — %;) for j = 1,...,d and k = 1,...,n, where a is the j*” component of x* and 77 
is the average of the j*” component of x’ for 7 = 1,...,n. Notice 


(40,43) = (Ldap = = (04 —H) = 


i=1 i=1 


het € = [eis ¢o5-0.64|, = |i, Ta," 2 3 Tal and y= So ui Then our approximation 
Gar 
more specifically takes the form 


KE x)= boi hatte —@;) =cotc:(x—X). 
j=l 
This form highlights that K (Fx) is a constant approximation with adjustments that 
increase in magnitude as x moves away from the center of mass. The coefficients c 
determine the orientation or tilt of the plane Kk. 
The goal of least squares is to minimize (K(F,x) — y, K(F,x) — y)q,, which is 
achieved at the solution to the following system of equations using the Gramian of 


the ; for 7 = 0,...,d. 


1 0 O sade 2% y 
C 

DO 8 ada See one |) cee 
Cc 

0 (odie Onde -- Ondae |]. [=| Wr dde 
Cm 

0 (ga, P1) p (ga, 2) p ie Oa: da) p (y, da) p 
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It is clear that co = y. Let us denote the real symmetric matrix A and vector b as 


follows 
(¢1, 91) p —_ (b1, ba) p 
A= : 7. : 
(bas 1) p as (a, Pa) p 
and 
(y; 1) p 
b i= (y, 2) p 
(y, Qa) p 


Then we solve Ac = b for c. As long as A is nonsingular, we can choose from several 
methods to solve this type of system. Caution must be taken when |A| is near 0, as 
our approximation blows up. If the measure was Lebesgue, there would be no issues 
with the matrix A. Missing data could cause the points in Ag to be concentrated in 
a small region away from v. Also, the data could lie on a lower dimensional space 
causing |A| to be close to zero. The truncated SVD method described in Section 
1.3 allows us to process the system Ac = b in all cases and model the appropriate 


geometry. The plane of best fit through the points F' is 
K(F,x) =Yyrt+e-(x—X). 


Consider the case when the points F’ are sampled from a power line (Figure 2.1). 
The system Ac = b is singular and the data in A has two singular values near or 
equal to zero. We can still solve the truncated system Ac = b for c. The rotation 
matrices U and V have been truncated to U; (3 x 1) and V,; (3 x 1), so the rotation 
of the solution c is limited. Thus, the dimensionality of the points F’ only affects the 
tilt of the plane K. 

In addition to the issue occurring when the determinant of A is close to 0, we also 


have a problem when (x — X) is large. The further away the x is from the center 
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Figure 2.1 A triangle where the points F' lie in a low dimensional space. 


of mass of F’, the less reliable our approximation becomes. Therefore a subdivision 
criterion limiting the distance between the center of mass of F' and the center of mass 
of Ag is added when creating the partition. 

Now we find the radius value rz at the intersection of Vv; and K for i=0,1,2. Let 
v = (0,¢) be any ray emanating from the origin. Then using the conversion formulas 


(2.3) and (2.4) we can write 


re cos(¢) cos(A) 
V =| recos(¢) sin(@) | for any re > 0. 
resin(?) 
Plugging these values in for AK (F’,x) we have the following equation to solve for re: 


resin(@) = y + c1(recos(¢) cos(@) — Fz) + c2(re cos(¢) sin(@) — XZ). (2.6) 


Then the intersection of the plane of best fit K and the v is 


Y — C1 XL — Co%Q 


sin(@) — c, cos(¢) cos(@) — c2 cos(@) sin(A) ” 


re= 


(2.7) 


The value rz is the local value of v for a triangle Ag. Hence we write r;(v, As). 
One potential issue is that A and v could be close to parallel. Then re would 
lie outside the range of the r-values in F’. The algorithm in Learning Theory [2] 


suggests that we clamp the value ry between ry := max{r : (0,¢,r) € F} and 
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Tm ‘= min{r : (6,¢,r) € F}. To clarify the meaning of clamp, when r; < r,, we set 


re =Tm. Similarly, if re > ry, then re = ry. 


Global Approximation and Discontinuities 


In this section we combine all of the local r-values for a vertex to a single r-value 
known as the global r-value. Each vertex may be contained in several different tri- 
angles. Let N(v) be the set of triangles in P, that contain the vertex v. Then we 
have the local values {r¢(v, As) }agen(v)- Note that in our construction of P4, N(v) 
contains at most eight triangles. 

For continuous regions of the partition we need to assign a single value to each 
vertex. This value is known as the global r-value of the vertex v and is denoted by 
r(v). We calculate r(v) as a weighted average of the local r;(v, Ag), where the weights 
are the sum of barycentric coordinates corresponding to v. Recall that A,(As, p) is 
the barycentric coordinate of the point p in the triangle Ag corresponding to the 
vertex v. Define the weight of the triangle Ag with respect to the vertex v to be 

W (v, Ag) := S- Ay(Ag, p) . 
pcAs 


Then the global r-value is 
S- ro(v, As)W (v, As) 
- AsEN(v) 


S- Wy, As) 


AséEN(v) 


r(v) 


In our particular application the approximation is not continuous everywhere, so we 
need to relax the continuity requirement. The global r-value is a weighted sum and 
does not approximate jumps well. For LiDAR scans, discontinuities occur when the 
distance from the scanner changes significantly over a small change in angle. Often 
this results from two different objects being sensed. Other occasions such as near 
the horizon are discussed in the remarks found in Section 2.3. If the discontinuity 


occurred within a triangle, then the variance of the r-values in the triangle would 
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Vx 


Figure 2.2. Two triangles A; and A, with shared vertex v*. The vertex is assigned 
two values. 


be large forcing the triangle to split based on the subdivision criteria. However, 
it is possible that the discontinuity is near an edge or vertex of the partition. A 
simplification of this situation would be two triangles A; and A, sharing a vertex v* 
and the r-values of Az are much larger than the r-values in A; (Figure 2.2). The 
solution is not just to subdivide further, because the same problem would just occur 
on a finer scale. Furthermore, both triangles could have small variance with respect 
to r, dense distributions, and have moderate changes in global r-values within the 
respective triangle, causing the triangle not to be marked for subdivision. Therefore 
our solution adjusts the vertex r-value, by allowing the vertex to have two different 
r-values. 

We determine the different values of r for a vertex v by clustering the triangles in 
N(v). For the simplification above, it suffices to separate N(v*) into Ni(v*) = {A} 
and No(v*) = {Ay}. The global r-value for v* is a weighted sum of rp(v*, A;) and 
re(v*, Ao). Now v* has one r-value r¢(v*, A;) for A; and another r-value rp(v*, Az) 
for Ay. It is possible to separate N into at most eight sets, as eight is the maximum 
number of triangles possible to contain a single vertex. It is a rare occurrence to 
need three or more sets to distinguish the objects. Therefore, we choose to separate 
N(v) into two sets Ni(v) and No(v) when v is along a discontinuity. If one wishes to 
separate into more sets, then the procedure is easily extrapolated. Once we determine 
the sets, then r;(v, N,) and ro(v, N2) are calculated by weighted averages just as we 


did before. 
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To determine the sets N, and Nz we march through the triangles in N. If re(v, As) 
is less than the global value r(v), then As € Nj. If re(v, As) is greater than r(v), 
then Ag € No. If re(v, As) = r(v), then Ag is in both N, and No. We only use 
additional values for r when it significantly decreases the variance around v. More 
specifically, we define the variance of N as the sum of the variances of the triangles 


in N. We decide to use these distinct r-values for v if 
Var(N) — Var(N,) — Var(N2) (2.8) 


is greater than some threshold. 


Algorithm Options and Remarks 


In this section we suggest some possible variations to the algorithm and make some 
remarks. We first discuss modifications to the subdivision criteria that produce bet- 
ter adaptive partitions, at the expense of computational efficiency. Another option 
discussed here is to vary the choice of the sets used to approximate the value at a 
vertex. Finally we give some remarks about special circumstances. 

The adaptive partition P, received is dependent upon the point cloud and the 
subdivision criteria implemented. There is a wide range of potential criteria. Here 
we outline some different options for the subdivision criteria that will serve the same 
purposes. In the description of the main algorithm, a triangle Ag is marked for 
subdivision if either Var,(Ag) or d(Ag) is larger than some respective threshold. 


Recall that Var,(Ag) can be expressed by the aggregate quantities of As, namely 


[i tdpa, | rdpa, and J. IrPapaw. 
As As As 


Moreover d(Asg) is expressed in terms of 


[spa | Pdpa. and |) 0a. 


40 


and the (0,@) coordinates of the vertices of the triangle. One desire is for P,4 to 
have triangles with reasonable distance between the center of mass of the triangle 


and the center of mass of the points in the triangle. An alternative would be to 


measure the distance of the centers of mass in R?. The issue here is that the 2, 


2, and y coordinates for the triangle vertices are not necessarily calculated. This 


requires that the r-values of the vertices have been calculated and converted into 


R? coordinates. Another possible subdivision criteria is to mark triangles with large 


changes in the r-values of the vertices. We do not include this in the main algorithm 
because it also requires the additional calculation of the r-values of the vertices at 
each subdivision level. The benefit of this criterion is that it reduces the number of 
long skinny triangles received along the line-of-sight of the sensor. In these variations 
the increased quality is counteracted by the decreased efficiency of the calculations on 
each level of subdivision. The criteria used in the main algorithm is faster, because 
only the angles are needed and they are known for all the vertices at the time of 
subdivision. A similar situation arises with the subdivision criteria requiring the 
variance to be small. The main algorithm uses the variance of the r-values of the 


points in the triangle for the criteria, but one could use the Cartesian coordinates. 


Using just the variance of y or any single dimension in R® as the criteria is not helpful 


because within one triangle any dimension could have a large variance. We resolve 


this by using the smallest eigenvalue of the Cartesian coordinate covariance matrix 


as a value describing the variance in R*. The additional calculation here is in finding 


the eigenvalue. 

When calculating the vertex r-value of a vertex v, a plane of best fit is found over 
some region containing v. Moreover, the main algorithm calculates a local r,-value 
over the triangles containing v and then a global value as a weighted average of the 
local values. An alternative is to fit the plane to a different neighborhood around v. 


Let U denote the neighborhood around v. Letting U equal X is too general and we 
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would lose the local topology around v. One idea: in the regions where the surface 
is continuous, we let U = N(v) and calculate the plane of best fit over U. In this 
case there would be no need for separate local and global calculations of r. For the 
discontinuous regions the set U could be clustered into U; and Uz as discussed in 
2.3. Then fit the plane to each U, and U2 and assign v the two resulting r-values. 
The downside to U = N(v) is that our approximation would be poor when v is 
topologically similar to a local maximum or minimum. Other options are possible for 
the choice of U and would affect the r-values and quantities computations accordingly. 

Now we discuss some issues that may appear in our resulting surface approxima- 
tion. The first is referred to as the horizon issue. A horizon is characterized by a 
flat surface nearly parallel with the sensor that extends far away from the sensor. A 
significant change in the radius may occur over a small change in the angle, forcing 
the triangles near a horizon to split because Var,(As) would be large. There are 
too few points over the large gap to conclude with any certainty that the surface is 
continuous and flat. Thus, we choose not to modify this criterion. Another issue is 
that the approximation of the vertex r-value is most accurate and stable at the center 
of mass of the neighborhood N(v). Thus we have the option to add a post processing 
step of moving interior vertices to the center of mass of the points contained in N(v). 
For clustered vertices, we split v into two vertices. One copy moves to the center of 


mass of Ni(v) and the other vertex copy moves to the center of mass of No(v). 


2.4 PROGRESSIVE TRANSMISSION OF AN APPROXIMATION 


We shift now to the problem of transmission of data. Not only do we want a quality 
surface approximation, but we also want an efficient encoding method. The multires- 
olution structure inherent in the newest vertex bisection method allows for varying 
levels of detail. The point cloud is considered to be very large and direct transmission 


of the entire point cloud is not reasonable. An example of the need for compression 
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is the case of approximating an unknown terrain. A robot or an unmanned vehi- 
cle may collect data, then transmit the point cloud information back to base. This 
transmission may progressively transmit the data. At first the robot may send very 
few bits, and then sends more as the user requests. Multiresolution method provides 
the framework for our encoding, and we introduce nonlinearity in the multiresolu- 
tion process itself. In our representation we give higher priority to certain geometric 
features (e.g. local extrema) by switching them with less significant elements in the 
multiresolution tree. We start by explaining our tree structure and then describe the 


switching process. 


Vertex Tree Structure 


Our encoding will be based on a vertex tree instead of a triangle tree structure. The 
idea is to build a tree to store information for the vertices. Let V be a tree where each 
node is a vertex. The initial partition Po has two triangles which form a diamond. 
Force-splitting the diamond results in the addition of a single new vertex v*. Then 
we let the root of V be v*. Splitting all of the edges opposite v* results in four new 
vertices, namely v1, V2, V3, and v4. Define the children of v* to be vj, V2, v3, and v4. 
Then we recursively define the children of a vertex v to be the new vertices added as a 
result of splitting the edges opposite v. Triangles of the partition are right-iscocleses, 
so the smallest angle is 7/4. This implies that no vertex can be contained in more 
than eight triangles. Therefore every internal node of the tree V has no more than 
four children. 

The tree V is very repetitive. Any vertex internal to the partition will appear 
twice, because of the force-splitting procedure. Whenever the partition is conformal, 
we have no need for this repetition. However, in the spherical approximation, some 
vertices have two distinct values and would make use of the duplication of vertices 


in V. So we define V’ to be the tree of vertices with root v* and when a new vertex 
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Vv; V, 


Figure 2.3. The vertex v* is created when the black diagonal edge E is subdivided. 
The red edges result from the subdivision of EF. The vertices v and vj, lie on EF, 
and the vertices v; and v2 are the vertices opposite E. 


is created, it is appended to the tree as the child of the vertex that initiated the 
subdivision. Thus V’ is a pruned subtree of V with no duplicate vertices. Whenever 
the subdivision is uniform, the tree V’ will be extremely lopsided with the majority 
of vertices on the left side. For uniform subdivisions, the tree structures for V,V’ 
and T are fixed base only on the dimensions of the initial square X. Therefore, our 


encoding of the surface will be based on these trees. 


Correction Value 


The location values of the vertices are automatically encoded within the structure of 
the vertex tree. Now we need to encode the y (or r) values of the vertices. Better 
than simply encoding the whole value, we predict a value of the vertex based on 
the four vertices of older generations (coarser level) around it. Then we only encode 
the difference, which we call the correction value. The correction bits are the binary 
digits of the correction value. We label these older generation vertices vi, v5, vi, and 
v2 as seen in Figure 2.3. 

Let C(v*) be the correction value of a vertex v*, and y(v*) be the predicted 
value of v*. Then C(v*) = y(v*) — 9(v*). The prediction is based on the y-values 


of the vertices vj, v2,v, and v5. The edge subdivided in the creation of v* is the 
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/ / 


edge from vi to v5. We let Maz = max {(y(v1), y(v2), y(v), y(v5)} and Min = 


min {(y(v1), y(v2), ¥(v1), ¥(v2) }- 


We have the following three variants of the prediction formula: 


Line Average Prediction 


av) = | 5(u(v) + u(v4))| 5 (2.9) 
Total Average Prediction 

av") = | F(u(vr) + ulve) + (v4) +y(vy))] 5 (2.10) 
Max and Min Average Prediction 

av") := F (Max + Min)| (2.11) 


Vertices on the boundary will not have the vertex v2 in the calculation of the predicted 
value. To modify the prediction formula for boundary vertices we remove v2 from all 
of the formulas and from the calculation of the Max and Min. Then in Equation 
(2.10) we multiply by 1/3 instead of 1/2. The formula for calculating the prediction 
value is not limited to the three we described here. These were chosen as simple 
predictions that are easily computed. For example, one could choose to use something 
similar to a plane of best fit through the vertices. Equations (2.9) and (2.10) are 


motivated by the L? error, and (2.11) is motivated by the L® error. 


Switches to Preserve Local Extrema 


In the standard multiresolution process, each stage of the progressive transmission 
would be the information corresponding to some subtree of V. Then the vertices or 
triangles not contained in the subtree would be approximated based on vertices in 
the subtree. One issue is that an important vertex may be present at a finer level not 
present in the subtree. To resolve this we introduce nonlinearity to the multiresolution 


process by switching the y-values of finer level important vertices with the y-values 


45 


of coarser vertices. For our purposes local extrema are considered important, as we 
aim to preform well in the Hausdorff error. This results in the coarsest level (Po) 
containing the global maximum and minimum y-values. As we refine, the previously 
switched vertices will be unswitched. 

Here we elaborate on the term switch. Vertices consist of a location and a y-value; 
however, in V vertices behave more like placeholders of the correction bits. All of the 
encoded information for v is stored in its spot in V. Switching the y-value of two 
vertices v, and vy is to switch the locations in V. So the correction bits for v, will be 
stored in the location of v, in V, and vice-versa. Furthermore the switch only occurs 
between vertices along the same line. One reason we do not switch with the other 
vertices, is that the generation gap between the vertices on the line is one where the 
gap is larger for the other vertices. Another reason is that the direction of the switch 
can be encoded with only one bit when there are only two choices. Let v4, be the 
vertex with the larger y-value along the same edge as v*, and vi, be the smaller of 
the two. 

We assume we already have a uniform partition to the finest level with the y- 
values of the vertices either given or approximated. Starting from the finest level we 
merge all of the triangles level by level until the initial partition is reached. Since 
our partitions are conformal we merge diamonds and boundary triangles instead of 
individual triangles. Each merge reduces the vertices in V by one and the triangles in 
T by one or two. The y-values of vertices are switched only when the vertex to remove 
is a local extrema. Recall from the prediction formulas the vertices v*, v1, V2, V1, V5, 
and the values Maz and Min. Now v* is the vertex to be removed instead of added. 
If y(v*) € [Min, Maz], then v* is not considered a local extrema. In which case, the 
vertex values will not be switched, and everything proceeds as normal. However, if 
y(v*) is outside the range of Maz and Min then v* is considered a local extrema and 


we switch the y-values. More specifically, when y(v*) > Maz the values y(v*) and 
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y(v,) are switched. If y(v*) < Min then the values y(v*) and y(v/,) are switched. 


Encoding and Decoding 


Progressive transmission of our surface approximation requires the following informa- 
tion to be sent: the triangle tree structure 7’, the vertex tree structure V, switches 
and correction values. Encoding the trees T’ and V only requires the dimensions of 
X and the finest resolution. The finest resolution is the length of the smallest edge 
in the finest partition, which is typically one. Therefore the header includes the di- 
mensions of X and the finest resolution. Other information included in the header 
is the four y-values of the initial partition and the index of the maximum correction 
bit. Note that the y-values of the initial partition contain the global maximum and 
minimum due to the switches. The correction values are sent bit by bit, so the index 
of the maximum bit allows us to know the magnitude of the first group of correction 
bits sent. We estimate the number of bits in the header to be about 84. 

The number of bits in the maximum correction value determines how many trans- 
missions will be sent. We refer to a transmission as a wave. The ith wave sends the 
information vertices with the 7th correction bit nonzero and all of their ancestors. The 
initial information for a vertex in V is encoded in either two or four bits. Consider 
the initial information for v. The first bit indicates whether or not v was switched, 
where 0 implies not switched and 1 means switched. If the first bit is one, then the 
second indicates which vertex v switched with. Recall that v could only be switched 
along the edge, so 0 and 1 are sufficient to indicate the vertex. On the other hand if 
the first bit is 0, then the second bit indicates whether or not we sent a correction 
bit for this vertex. If the first two bits are 00 this means that no information about 
v nor any vertex below v in V is sent, and is called a zero-tree. Signifying a zero-tree 
tremendously reduces the number bits sent on each wave, especially with the early 


waves. If not a zero-tree, then the third and fourth bits are sent. The third bit 


AT 


Figure 2.4 A SICK Commercial LiDAR scanner. 


indicates the sign of the correction value and the fourth bit is the ith correction bit. 
Note that the initial bits for vertices far down the vertex tree V are often sent in a 
later wave. For v the waves after the one with the initial bits will only contain the 


corresponding correction bit and nothing else. 


2.5 NUMERICAL RESULTS 


Virginia Tech Data 


Andrew Kurdila, Department of Mechanical Engineering & Unmanned Systems, Vir- 
ginia Tech and his team provided the three point clouds analyzed in this chapter. Each 
point cloud was created from two SICK Commercial LiDAR scanners (www.sick.com) 
mounted on a vertically rotating motor, which the team registered and combined to 
form the data sets (Figure 2.4). The SICK scanners are used by ARO Multidisci- 
plinary University Research Initiative (MURI) and USAT. The scanners are consid- 
ered stationary. The point cloud VT was created from two stationary scans of the 
courtyard of Hancock Hall, Virginia Tech in connection with MURI Topic # 28, Dy- 
namic Modeling of 3D Urban Terrain, which one can access at the following website: 
http: //imi.cas.sc.edu/MURIwebsite/aro-broad-agency-announcement. 

The 81,150 points making up the point cloud VT are provided in spherical co- 
ordinates (longitude, latitude, distance). The color scheme for the spherical point 


cloud and approximations range from green to red as the latitude increases and blue 
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is added as the distance from the sensor increases. This scheme results in an intuitive 
visualization of the spherical point cloud. Standard coloring schemes color according 
to the r-value, z-value, or the intensity value given in the LiDAR scan. The detailed 


RGB color is 


(S42P 1 StaP r 


1 nm  °max{r} 
for the coordinate (0, ¢,r). 

Figures 2.5 and 2.9 are different views of the original VT point cloud consisting 
of 81,150 points, and are both distortions between Cartesian space and spherical 
space. The other figures are of the approximation generated by the algorithm in 
this chapter. Figures 2.10 and 2.11 are of particular interest as they highlight the 
phenomenon that occurs along the line-of-sight from the sensor. The objects sensed 
in VT include three trees, a building and a hallway. The algorithm was executed 


with the following parameters from Section 2.3: 
e Splitting Criteria: Triangles with Var,(As) > 0.8 or d(Ag) > 0.5 are subdivided. 


e A vertex will have two distinct r-values when 


Var(N) — Var(Ni) — Var(N2) > 3. 


Compression Results 


In this section we present the results of our compression on a 2049x2049 STM (stan- 
dard terrain map) file of the Grand Canyon. The maximum height value is 1709 
meters and the minimum is 1336 meters. Hausdorff error is reported in meters. We 
compare our results in Tables 2.2, 2.3 and 2.4 with that of JPEG 2000 using the 
Hausdorff error, which is found in Table 2.1. Naively we send our bits directly and do 
not compress the bitstream itself. The team who developed JPEG 2000 invested time 


and effort on the compression of their bitstream. Therefore, it is not appropriate to 
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Figure 2.5 Slightly distorted view of Figure 2.6 Surface generated from 


the VT point cloud between spherical the algorithm. From left to right a 
and Cartesian coordinates. (81,150 hallway, three trees, and a building 


pts) are distinguishable. (13,781 triangles) 


rs 
Figure 2.7 Improved the Figure 2.8 Moved the interior 
approximation by adding the optional vertices to their respective centers of 
subdivision criterion of maximal mass to improve the accuracy of the 
radial distance 2. (17,566 triangles) approximation. (17,566 triangles) 


directly compare our method and JPEG 2000. We provide the following as ballpark 
measurements of performance. Our results show that for extreme compression and 
full recovery our encoding outperforms JPEG 2000. For intermediate Hausdorff er- 
rors between 3 and 10, JPEG 2000 is slightly better. JPEG 2000 was not particularly 
geared for Hausdorff error, and we are not surprised when it does not perform well 


in Hausdorff error for the extremes. 


Table 2.1 Wavelet Results (JPEG 2000) 


Bits Original Bits Compression Hausdorff Error 
1 || 2576 67174416 .00383% 52 
2 || 14024 67174416 .02088% 8 
3 || 61512 67174416 .09157% 3 
4 || 80912 67174416 .12045% 3 
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Figure 2.9 Notice the occluded Figure 2.10 Long skinny triangles 
points behind the large tree along the occasionally appear along the 
line-of-sight. line-of-sight. 


Figure 2.11 ‘Triangles with radial 


distance larger than 2 are further Figure 2.12 Moved the interior 
subdivided, so there are much fewer vertices to their respective centers of 
skinny triangles along the mass to improve the accuracy of the 
line-of-sight. approximation. 


Table 2.2 Line Average Prediction 


Triangulation Bits Original Bits Compression Hausdorff Error 
1 390 67174416 0.00058% 18 
2 2834 67174416 0.00422% 15 
3 18008 67174416 0.02681% 12 
4 75689 67174416 0.11268% fi 
5 309623 67174416 0.46092% 4 
6 1136636 67174416 1.69207% 2 
7 3358479 67174416 4.99964% 1 
8 7648647 67174416 11.38625% 1 
9 14620231 67174416 21.76458% 1 
10 26335126 67174416 39.20410% 0 
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Table 2.3 Total Average Prediction 


Triangulation Bits Original Bits Compression Hausdorff Error 
1 150 67174416 .00022% 27 
2 1700 67174416 .00253% 16 
5 12168 67174416 .01811% 11 
4 53434 67174416 .07955% 7 
5 221580 67174416 .32986% 4 
6 862531 67174416 1.28402% 2 
7 2775789 67174416 4,13221% 2 
8 6801281 67174416 10.12481% 1 
9 13858297 67174416 20.63032% 1 
10 25767964 67174416 38.35979% 0 


Table 2.4 Max and Min Average Prediction 


Triangulation Bits Original Bits Compression Hausdorff Error 
1 242 67174416 .00036% 30 
2 1637 67174416 .00244% 17 
3 13473 67174416 .02006% 10 
4 62768 67174416 .09344% 6 
5 260524 67174416 .38783% 4 
6 1013978 67174416 1.50947% 2 
7 3163858 67174416 4.70992% 1 
8 7454856 67174416 11.09776% 1 
9 14528905 67174416 21.62863% 1 
10 26218411 67174416 39.03035% 0 
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CHAPTER 3 


DISTRIBUTION-DEPENDENT SUBDIVISION SCHEMES 


For a point cloud D, the points are divided into subdomains (cells) based on the 
partition. To improve the computational cost, we wish to only have the aggregate 
quantities Q representing the points in each subdomain. There are two types of quan- 
tities in Q: distribution-dependent Qx and value-dependent Qy. The distribution- 
dependent quantities depend only on the location and distribution of the points in D. 
The value-dependent quantities involve the y-value of the points in D. When a cell of 
the partition is subdivided, Q of the children needs to be calculated based solely on 
the Q of the parent and the neighbors of the parent. The general idea of subdivision 
schemes is to find the function values of finer vertices using the values of vertices from 
coarser levels. Unique to our algorithm is that our quantities are both distribution- 
dependent and value-dependent. The method described in this chapter addresses the 
issues of irregular sampling and gaps in point clouds. Traditional methods operate 
under the assumption that a single cell from subdivision has Lebesgue measure, how- 
ever we wish to learn the measure within each cell as well as the function value. One 
particular benefit of learning the measure is that cells on the boundary are treated 
the same as cells with neighbors of zero measure. No special scheme is needed for the 
boundary cells. 

We first approximate the distribution-dependent quantities Qx then use them 
to construct a polynomial approximation of the value-dependent quantity Qy. The 
foundation for this approach is explained over the intervals (d = 1) in Section 3.1. 


Then in Section 3.3 we discuss the algorithm over triangles (d = 2). 
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Recall that D is a point cloud with distribution p, and marginal probability mea- 


sure px. Over an interval J € X, let 


Qx (UU, px) = (MoU, px), Mi, px), Mo, px)) (3.1) 
= (/ ldpx. [ edox, | 2*dpx) 
Qv(L,px) = My(L,px) = f vdpx. (3.2) 


For a triangle A € X, let 


Qx(A, px) = (Moo(A, px), Mio(A, px), Moi(A, px)) (3.3) 
— (/, ldpx, f wdpx, f wdpx) 
Qy(A, px) = My(A, ex) = f dex. (3.4) 


Optionally, one could use the Bernstein polynomials instead of the traditional 


moments and let 


Qx(L.ex) = (f Bile)dox, | Bie)dox, | B3(e)dpx) (35) 


and 


Qx(A,px) =(f, Bholdex, | Binodex, f Binder). (8.6) 


Note that the quantities based on the Bernstein polynomials are also additive. The 
use of the Bernstein polynomials increases the numerical stability. Furthermore, 
the Bernstein polynomials are simply a change of bases from the standard poly- 
nomial bases. So, one may transform the quantities between the two bases as de- 
sired. To simplify the explanations we use the moments in the standard bases for the 
distribution-dependent quantities in the descriptions and calculations of Sections 3.1 


- 3.3. 
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Since px is unknown, we approximate the moments by summing over the points 


in each cell. For the points {(x;, y;)}"_, in the interval I we have 


and for the points {(x, 7, y’)}™, in the triangle A we have 


Qx(A, ox) ones So) (3.9) 


3.1 APPROXIMATION OF Q OVER INTERVALS 


Approximation of the Distribution-Dependent Quantities 


Standard subdivision schemes reproduce the polynomials on each subdivision step. 
As an analogy to polynomial reproduction, at each step we want to reproduce the 
moments for certain types of measures. Given the moments of an unknown measure 
p over some region R, we wish to find a known measure p such that Qy(R,p) = 
Qx(R, 4), and we say that ys represents p. Then we approximate the distribution- 
dependent quantities of the children with respect to p by the direct calculation of 
the quantities with respect to the known measure ps. The representative measures we 
consider are parameter-dependent measures so they can easily be determined. 


First consider probability measures jz of the form dy = w(x)dx with: 
e Polynomial Weight Measures: 
w(r)=a+ Bret yx? 


The parameters of this type of measure are (a, 3,7). Since yw is a probability 


measure, there is a constraint equation which defines a in terms of 3 and y¥. 
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e Subdomain Measures: 


1 — X{a,8] 
Opa. ate) = Coe) 


The parameters of this type of measure are (a,b). These measures are Lebesgue 
measures of either [a,b] or [—1, a] U[b, 1] for -1 < a < b < 1 and single point 


measures when the Lebesgue measure is zero. 


Once we establish the schemes over these types of probability measures, then the 
results are extended by introducing a constant C’ such that jy is now of the form 
du = Cw(x)dx. The main idea of the method does not change with the addition 
of the constant, there is just one more parameter to calculate before we know the 
measure. In both cases above, the additional calculation is simple. Our final scheme 
blends the two schemes together because each scheme on its own does not make for 
a good method, but together they perform well. Polynomial weight measures do not 
have the range to be able to represent all possible measures. On the other hand, we 
show that measures can always be represented by a subdomain measure, but they 
introduce a gap for any measure not exactly Lebesgue on the entire interval. This 
is unlikely to be realistic for measures close to Lebesgue. Thus, we combine the two 
schemes into one blended scheme. 

We describe our schemes over the interval [-1,1] as there is an affine transformation 
to and from a general interval and the interval [-1,1]. Thus, the parent interval is 
[—1,1] and is divided into its left child [—1,0] and its right child [0,1]. Here we 
describe the details of the transformation. Any element from a general interval |r, s] 


can be transformed to the standard interval [—1, 1] by the affine transformation 


2 
re) = eee 


S—Tr r—s 


whose inverse is 


Using the transformation T we also transform Qx (|r, s]) to Qx({-1, l)). 


Mo([-1,1]) = : e ~Mollr, s]) (3.11) 
Ma((-A, 1) = Mlle sl) — 2 aol, (3.12) 
: a A eee RE) 


Mo( —1,1 ) = eae s]) _ (s—r)3 


(r — s)3 
If Bernstein polynomials are used for the moments instead of 1,7, and x”, then the 
quantities would be affine invariant and this transformation would not be needed. 


We assume during the development of our algorithm that the interval is [—1, 1], then 


using the transformations above, arrive at the quantities over the general interval. 


Description of the Geometry of the Moments 


Here we establish the setup for our geometric reasonings. Let fz be a probabil- 
ity measure, then M)({-1,1],u) = 1. Define u(u) := M,([-1,1], 4) and v(m) := 


Mo2([-1, 1], «). Thus the quantities Qx([-1, 1], u) = (1, u(y), v(u)) and each measure 


4. corresponds to a vector q(j1) := (u(j),v()) € R?. We refer to this coordinate 


system as the uv-plane. The children’s quantities can also be expressed in the uv- 
plane by restricting ys to [—1, 0] or [0,1]. Then Qx([—1, 0], “) = Qx([-1, 1], H]-1,0))- 
For the purpose of the uv-plane we ensure it is a probability measure by dividing 


Mo([-1, 0], 4). In fact, for any measure ps we define 


— M(t U4) = 1] 
MON (A, Ta) 8 = Nig(=A Tha) 


and again q(u) = ((u(), v(4)). The vector qg is not unique to the measure py. In fact, 
there are infinitely many probability measures such that the corresponding q is the 
same vector. From the parent’s vector in the uwv-plane, our algorithm approximates 
the children’s uv-vector. 

Now we explore the geometric properties of the wv-plane and the coordinates of 


q(u). The Dirac delta function can be considered a probability measure concentrated 
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j i H 7 Figure 3.2. With respect 
to moments based on the 


Figure 3.1 The curve g is the red curve in the Bernstein polynomials, Q 
uv-plane. 2 is the convex hull of g. The yellow is the convex hull of the 
dash line indicates the convex hull of the left and red curve. Each axis 

the right of g. The vectors corresponding to the corresponds to the 
children will lie between g and the yellow dashed respective Bernstein 

line. moment. 


at a single point x in the interval. We refer to these delta functions as single point 
measures 6, for some x € [—1,1]. Then we have that Qx({-1, 1], 6.) = (1, 2,27), so 
u(d,) = x and v(6,) = x”. Let the curve g be defined by g(u) = u? for all u € [—1, 1] 
in the uv-plane. Then q(6z) = (u(d.), g(u(d.)) for all single-point measures 6,. We 
define 2 to be the convex hull of g. Then the vector q for any probability measure 
is contained in 2. Loosely speaking any probability measure can be considered as a 
convex combination of the corresponding single-point measures. 

The main idea on the interval [—1,1] is to use the distribution-dependent quan- 
tities Qx([-1,1],) and the geometric behavior in the wv-plane to approximate 
Qx({-1, 0], e) and Qx([-1, 0], p). The distribution-dependent quantities are very lo- 
cal quantities, so Qx([-1, 1], ) are the only quantities used to find Qx({—1, 0], p) 
and @x([0,1],). To do this a one-to-one correspondence is established between the 
measure parameters and the appropriate region in the wv-plane. However, an explicit 
formula may not be known, in which case we develop a scheme for finding (u,v) of 


the children. 
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Polynomial Weight Measures 


In this section we wish to reproduce polynomial weight measures 7 of the form dn = 
Cw(x)dx where w(x) = a+ Bx + yx?. We first develop the method for probability 
measures, then calculate the constant C’. Considering 7 a probability measure of the 


form dn = w(x)dx produces two properties: 


e i (a+ Bx +72?) de =1 


-1 


e w(r) =a+8x+ ya? > 0 for all x € [-1, 1]. 


Thus 
a = 1/2 — 1/37. (3.14) 


Ensuring that the measure is nonnegative produces more restrictions on a, 3, and 
y. Thus w(x) > 0 for all « € [-1,1]. More specifically w(0) = a > 0. Then using 
the equation (3.14) we have 1/2 —1/3y > 0 => y > 3/2. We also have that the 
endpoints must be nonnegative, soa+6+y7y2> 0anda—8+y7 > 0. Adding 
the two together we have a+ 7 > 0, which implies that y > —3/4. We also get 
that [8] < a+y7 = 1/2+ 2/3y. To summarize, we have —3/4 < y < 3/2 and 
|B] < 1/2 + 2/3y. 

When —3/4 < 7 < 0, we have that w is either a concave down parabola or a 
line, so we only need to satisfy the endpoint criteria. Then the current restriction 
of 3 is sufficient. The restriction on ( is not sufficient when the w is a concave up 
parabola with the minimum inside of the interval [—1, 1]. Now consider y positive, 


we have that w is a concave up parabola with a minimum at 7 = Plugging the 
Y 


constraint on a and the location of the minimum into w, we have the minimum value 
6y — 477 — 36? _ 
Say S8ee gee 
127 27 


4 
or equal to 0. So 6y—47? — 36? > 0, which implies |6] < ,/2y — mile In the standard 


[—1, 1], then the minimum value needs to be greater than 
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conic form we have, 


4, 16 3\2 
3h 4 i (y- ) ahs (3.15) 


3 3 
which is an ellipse centered at y = Z and @ = 0 with major axis Z along y and minor 
3 
axis > along (. 
We now determine which values of y require the constraint (3.15). For 0 <7 < 


3/8, we have the following string of inequalities 


| 4 
1/24+2/3y > 4/27 — ane > 24. 


If |G| > 27 we only need the end point criterion |] < 1/2+2/3y as the minimum is not 
n [-1,1]. If |G| < 27, both criteria are already satisfied. So when —3/4 < 7 < 3/8, 
we have that |G] < 1/2+2/3y is sufficient to ensure w is nonnegative on [—1, 1]. The 


remaining values 3/8 < 7 < 3/2, give us the following string of inequalities 


4 
Qy > 1/24 2/37 > 4/27 - a 


Thus || is forced to be less than 27 by the end point criterion. Then we have that 


| 4 
\C| < 4/2y- an is the appropriate constraint on 8 given 3/8 < y < 3/2. Notice 


4 
that for y = 3/8 the value ,/2y — ar = 1/24 2/37. 


To summarize, there are two regions in the parameter space of 3 and y. The first 


region R, is given by 
a0 3 2 
—<vy< = d <=+ 7. 
(Se Cond? IBIS Gea 


The second region Rg is given by 


4 1 2 
eS7S5 and 364+ (4-9) <1. 
8 2 
In the uv-plane we have that 


Fait gee (3.16) 


Figure 3.3 The blue region denotes the region [in which the vectors 
corresponding to polynomial weights lie. The vector (0, 1/3) corresponds to 
Lebesgue measure over the entire interval [-1,1]. 


So the image of the region R, in the uv-plane is 


1 5 1 
— and |u| < -v—-. (3.17) 
5 2 2 


The image of the region R2 in the wv-plane is given by 


2 3 1542 Th 
ss. and 3u? + (>) (v- =) <1. 


We let [ be the region in the uv-plane corresponding to the transformations of both 
R, and Ry. Loosely speaking, I is the set of wv vectors for polynomial weights. Now 
consider the shape of the region I’. One region is defined by an ellipse centered at 
(0, x) with major axis ve along u and = along v. The other region is given by the 
interior of the absolute value function v = 2/5|u| + 1/5. The boundaries of these two 
regions intersect at (+5. =). In fact, the lines v = 2/5u+1/5 and v = —2/5u+1/5 
are the tangent lines to the ellipse at the points (+5. =). To show this we examine 


2°5 
the implicit derivative of the ellipse, which is 


15? 7\ du 
Out (o-—) =e 


dv —Au 12 dv 2 
This implies th = h I ea) have that — = +-. 
is implies that du 75v — 35 and at the points ( 9 =| we have tha a 5 


The correspondence between (3, y) and (u,v) is one-to-one and the inverse relation 


is 
45 15 


3 
PS5 and 9 er ae (3.18) 
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Now we formulate and prove the result from this section. 


Lemma 3.1. Let p be a measure on |-1,1]. If q(p) in T, then we construct a poly- 


nomial weight measure n of the form dn = Cw(x)dx such that Qx([-1,1],p) = 
Qx((-1, 1], n). 


Proof. The inversion formulas (3.18) give us the explicit values for G and y. From ¥ 


we use (3.14) to determine a. So far we have 


u(p) = i ah = I cw(x)dx and (3.19) 
1 2d 1 
u(p) = o Te = 2 xw(ax)dx (3.20) 


Thus we let C = f', ldp. Then Mo([—1, 1], 2) = C = Mo({-1, 1], n). Furthermore 
1 
My([-1, 1}, p) = u(p)Mo([-1,1],p) =C f ew(2)de = Mi((-1,1],) 
and finally 


Ma({-1,1],p) = v(p)Mo([-1.1],9) =C f° 2?w(x)de = Ma([-1, 1], 


Therefore, Qx([-1, 1], p) = Qx([-1,11,n). 


Given a measure p with q(p) in’, we determine the representative measure 7. The 
representative measure 7 is completely known at this point, and we simply calculate 


the quantities for the children. 


0 0 0 
Qx([{-1, 0], e) & (/ Cu(x)dz, | Caw(x)de, | Cx*u(x)dr) (3.21) 
= -1 = 

1 1 1 
Qx((0, 1], p) & (| Cuu(x)dz, | Cow(x)de, [ Cx°w(x)dr) (3:22) 

0 0 0 
If the measure 17 is precisely of second degree polynomial form, then our approximation 
is perfect. Because the acceptable region I’ is not equal to 2, some measures will not 


have a polynomial weight representation. Thus we need to develop an additional 


scheme. 
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Subdomain Measures 


In this section we construct a scheme that preserves subdomain measures of the forms 


dus = Cw(ax)dx where 


wot) = or w(r) = (1 =xieai) 


b-a 2—b+a 


Measures of the first form are referred to as interval measures and of the second form 
as complementary measures. As we did for polynomial weight measures, we first de- 
velop the method for probability measures, then calculate the constant C’. The main 
result from this section is that subdomain measures of these forms create an equiva- 
lence class for all measures. Two measures p and yu are in the same equivalence class 
if Qx({-1, 1], e) = Qx({-1, 1], w). In other words, given a measure p, we construct a 
measure py of the form above such that q(p) = q(z). We then calculate the children’s 
quantities over js as an approximation of the quantities over p. Once the parameters 
C’,a and b are determined, then we know the measure and thus know the measure 
of the children. As it turns out, finding a and 6 is challenging for these measures so 
we first prove that a one-to-one correspondence exists and then develop a scheme to 
directly find (u,v) of the children. 

First we consider the interval measures j1 of the form du = *“! on [—1,1]. The 
possible values for a and b are —1 < a < b < 1, which is a triangle in the ab- 
plane with vertices (—1, 1), (—1, —1) and (1,0). To clarify, when a = b the measure is 
concentrated entirely at a. In this case, w := dg. When a = —1 and b = 1, then pis the 
Lebesgue measure over the entire interval [—1, 1]. Let S = {(a,b): -l<a<b< 1} 
and let S° denote the interior of S. We let M = M(a,b) = (u(a,b), u(a,b)) denote 


the transformation from the coordinates (a,b) € S to the wv-plane by the equations 


b 
1 ade a+b 
u(a, b) = [edn = - ay (3.23) 
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Figure 3.4 On the left we have the set S in the ab-plane. On the right is the 
corresponding region A in the wv-plane. The colors indicate where the region and 
boundaries are mapped by M. The region corresponds to interval measures. 


and 


b 
2 
fed oe ader 
b-a 3 , 


v(asb) = I. xdu = (3.24) 


This region in the uwv-plane is referred to as A, and M: S > A. Note that M(a, b) 


corresponds precisely to q(j) when pu is of the form du = Meth dey, 


(1—xta,0}) ane 


Next we consider the complementary measures ut of the form dl = “*sae 


on [—1,1]. When a = 0, then df = 1/2dx. The only limiting case occurs when 


(a,b) € S approaches [—1, 1]. 


Lemma 3.2. As (a,b) € S approaches [—1,1], M"(a,b) is on the line v = 1 in the 


uv-plane. 


Proof. Consider the case a = —1 and 6 approaching 1. Then we have that oT 4, b) is 


1-2 1-063 
LS a1 b+)’ which goes to 1 as b approaches 1, and v = 3(<1-b +2)’ which 
also goes to 1 as b approaches 1. One the other hand, when b = 1 and a approaches 
2). 1 3 1 
—1, Ma, l)isu= ea which goes to —1 as a approaches —1, and v = any 


which goes to 1 as a approaches —1. So the coordinate u depends heavily on the way 


in which we approach (—1,1). Now for 7,¢ > 0, consider a = —1+.€ and b=1-—Te. 


T-1 
As € approaches 0, then wu approaches Tq and v approaches 1. Notice that u can be 
T 


any value in (—1,1) as 7 is an arbitrary nonnegative constant. Therefore, (a,b) € S 


approaching (—1,1) implies that M aCe b) is on the line v = 1. 
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Figure 3.5 On the left we have the set S in the ab-plane. On the right is the 
corresponding region Y in the uv-plane. The colors indicate where the region and 
boundaries are mapped by M®. The region corresponds to complementary 
measures. The limiting case (-1,1) is mapped to the line v = 1 and the hypotenuse 
a = b is mapped to (0, 1/3). 


Recall that S° = {(a,b) :1<a<6< 1} is the interior of S in the ab-plane. Let 
M° be the transformation from the vector (a,b) € S° to the wu-plane with respect to 


the complementary measures. Then MW Ca, b) is given by the equations 


a 1 
‘i adx + xdx 2 2 
_ ae b —_ «a =o 

u(a,b) = f edu ~" 9-b+a  22—b+a) 20) 

and ; 

2d d 
ite Pte rt foe ~ _ 2-B +08 (3.26) 
, i 2—b+a 3(2—b+a)— 


Let TY be the region given by (3.25) and (3.26) in the uv-plane. 
To better understand the regions A and Y we examine the transformations of the 
boundaries of S. Table 3.1 gives the M and M® transformations of the boundaries 


in S to the wv-plane. Define the boundary curve go(u) = 1/3(4u? — 2|u| +1). 
Lemma 3.3. A is bounded by the curves g(u) and go(u) in the uv-plane. T is bounded 
by v(u) = 1 and go(u). 


Proof. By definition of the regions A and Y, M(S) = A and M*(S°) = Y. Consider 
the transformations of the boundaries for S. 
M(a,1) = (1/2(a + 1),1/3(a* + a + 1)), which corresponds to go(u) for u € 


[0,1]. Then M(—1,b) = (1/2(b — 1),1/3(b? — b+ 1)), which corresponds to go(u) 
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Table 3.1 Boundary Transformations of S in ab-plane to the wv-plane 


S Bounds | Transformation M Image in wv-plane 
b= 1 yes = el ie Se u € [0,1] 
ne abd ere _ su? +a $1 é[-1,0 
— eae sa : u(u) =u? €[-1,1 
S Bounds | Transformation M® Image in wv-plane 
eal waist ot vu) = ett u € [-1,0] 
b+1 bP +b+1 4u? — 2u +1 
a=-1 oes ee a u € (0, 1] 
Ga) = 0) p=173 (0, 1/3) 


for u € {[—1,0]. Furthermore, lim M(a,b) = (a,a”), which is the curve g(u). The 
a, 

transformation M is continuous; therefore, M(S') = A is bounded by g(w) and go(u) 

in the wv-plane. 


Now consider the complement transformation M"(a,b) = (u,v) where 


for any (a,b) € S°. M C is not defined for any of the boundaries of S, so we will examine 
the limits. By Lemma 3.3 we have that (a,b) — [—1,1] gives us v = 1 in the wv- 
plane. Then lim M(a,b) = (1/2(a — 1), 1/3(a? — a +1), which corresponds to go(u) 
for u € [—1,0]. We also have the limit Jim, M(a, 6) = (1/2(6 + 1), 1/3(6? +6 +1), 
which corresponds to g2(u) for u € [0,1]. And finally lim M(a, b) = (0,1/3). Since 
M® is a continuous transformation on S®, we have that the region M CS °) is contained 


within the region bounded by gp and v = 1. Therefore, YT is the open region bounded 


by gz and v = 1. 


Lemma 3.4. Let p be a probability measure with q(p) in the interior of Q. Then 


there exists a measure yt of the form du = w(x)dx with 


= a 
Oe) : = _ 


_ Xtal op 


OG 5 
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Figure 3.6 The union of A and Y. Notice that together they cover the convex hull 
of g. 


such that q(p) = q(t). 


Proof. From Lemma 3.3 we see that AU YT is the interior of Q and that AQ Y = @. 
Essentially, we need to show that the transformations M and M® are invertible be- 
cause these transformations directly relate the representative measure parameters as 
a uv-vector. This further simplifies to showing that the Jacobian of the transforma- 
tions are nonzero on the domain. The Jacobian of the transformation M is 6(b — a), 
and only equals zero at a = b. Since the Jacobian is nonzero for (a,b) € S, we have 
that M is invertible. 
Now consider the transformation M", whose Jacobian is 


(b — a)((a — b)? — 12ab — 4) 
6(2 —b+a)3 


J= 


We want to show that J is not zero for (a,b) € S°, so we need only consider the 
function f(a,b) = (a — b)? — 12ab — 4. Note that f(—1,1) is zero, but (—1,1) is not 


in the set S°. We have that the gradient of f is given by 
V(f) = (3(a— 6)? — 126) i+ (-3(a —b)? — 12a) 3, 


and is zero at (-1,1). Thus f has a critical point at (-1,1) and f(—1,1) = 0. The 


second derivatives are 
fac = Te = 6a =) and fas = 6(b — a — 2) ; (3.27) 
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which are all defined and less than zero for (a,b) € S°. By the second derivative 
test, f(a, b) has a local maximum at (-1,1), and the second derivatives continue to be 


negative for all (a,b) € S°. Therefore, J is nonzero on S°, implying the transformation 


M° is invertible. Therefore, the transformations over A and Y are invertible. 


In our setup we are given a vector (u,v) and we want to find the corresponding 
(a,b). So we now know that such relation exists for (u,v) € 2. However, there is no 
explicit formula for the inverse of M and M°. Thus, we develop a scheme for finding 
the wv-coordinates of the children. First we develop the scheme for vectors wv on the 
boundary, then for uv on the interior of 1. 

Our main concern is not about the measures corresponding to the boundary, but 
we include the scheme here for completeness of our algorithm. Denote the wv-vectors 
of the children by gz and qr. The first boundary case is when (uw, u?) is the wv-vector 
corresponding to the measure. Then the representative measure is 6,,. If wu. € [—1, 0], 
then gz = (u,uj) and gr = (0,0). Similarly, if wu, € [0,1], then gz = (0,0) and 
dr = (u1,u7). The second boundary case occurs when the coordinates are (u2, 1) so 
v = 1, which corresponds to the limiting case. Then if wg € [—1,0], we have that 
di = (ua, 1) and gr = (0,0). Similarly, if ua € [0,1], then gz = (0,0) and gr = (ua, 1). 

The scheme for the interior of 2 requires that we further partition the sets A and 


Y. The additional boundary curves 


(u) —Alu|? + 4u?,/u? + 4ful — 12u? + 4]u|,/u? + 4|u| — 2 
93\U) = 
—3]u] + 3,/u? + 4|u] — 6 


and gi(u) = 4/3u? are needed. As it turns out, these curves distinguish the case 


ab > 0 from the case ab < 0. The curve g, divides A and the curve g3 divides TY. 
Essentially, we have one scheme for a and 6 in a single child and another scheme for 
a and b in different children. 

We define A; as the region bounded by g, gi, and ge including g, from —1/2 to 
1/2 and gz for u € [—1,—1/2]U[1/2,1]. Let A» be the region bounded by g; and 
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Table 3.2 Additional Transformations in the ab-plane to the wv-plane 


ab-plane | Transformation M Image in wv-plane 
b=0; |las=a/2 w= E13 v = 4/3u? u € [-1,0] 
a=0. | w=5/2 U0 /3 v =4/3u? u€ [0,1] 
ab-plane | Transformation MM Image in wv-plane 
b=0 Ga ie (u) é [0,1] 
— — — — U 
“a4 ° 3q-+6 ee 
0 J 2 (u) we [-1,0] 
a 3b — 6 a 


g2 including go from —1/2 to 1/2. Thus, A = A, UAg. Similarly, define TY, to be 
the open region bounded by go, and g3. Lastly, define Y2 as the region bounded by 
v(u) = 1,92 and g3 including the curve g3 from —1/2 to 1/2. So TY = TY, UYo. For 
A in the following lemmas, we divide S into the regions R; = {(a,b) € S': ab > 0} 
and Ry = {(a,b) € S: ab < O}. For Y, we have T; = {(a,b) € S° : ab > O} and 
To = {(a,6) € S°: ab < O}. 


Lemma 3.5. M(R,) = Ay, M(R2) = Ao, M*(T;) = YT; and M* (Tp) = To. 


Proof. By the definition of A and Y, we already have that M(S) = A and M°(S°) = 
Y. Consider the boundary a = 0 for all the regions. Then M(0,b) = (b/2, 67/3) 


which is gi(u) for u € [0,1], and 


Met0,b) = ( b? = 


2b—4’ 3b-6 
which is g3(w) for u € [-1,0]. When b = 0 we have that M(a,0) = (a/2,a?/3) which 


is gi(u) for u € [—1,0], and 


oe GD 
M®(a,0) = 
ae) (a 3) 


which is g3(u) for wu € [0,1]. Recall the boundary transformations in the previous 
lemma. 
Both M and M® are continuous transformations on their respective domains, so 


we have that the image of the subregions of the domains is contained in the map of 
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Figure 3.7 The top left is the set S in the ab-plane colored according to the 
transformation M and the top right is colored according to the transformation M C 
The bottom plot is the image in the uv-plane of the four regions A; (blue), Ag 
(pink), YT; (green), and Y2 (yellow). 


the boundaries. All that remains is to test one (a,b) transformation for each region. 
Then M(1/4, 1/2) = (3/8, 7/16) and 7/16 < gi(1/4) which implies that M(R,) = Aj. 
M(-1/2, 1/2) = (0,1/4) and 1/4 > gi(1/2) which implies that M(Rz) = Ao. 

Next M"(1/4,1/2) = (—3/56, 121/336) and 121/336 < g3(1/4) which implies 
that M°(T,) = Y1. M®(—1/2,1/2) = (0,7/12) and 7/12 > g3(0) which implies that 
M"(To) = To. 


Algorithm 


Here we describe the nonlinear algorithm we developed to find the wv-vectors of the 
children, given the wv-vector of the parent [-1,1]. The outline for the algorithm is as 


follows. Let p be any probability measure. We first determine in which region q(p) 
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lies. If g(p) is on the boundary of 2, then we follow the scheme discussed above. 
Thus q(p) is in one of the regions Ay, Az, 1, or Tz. The region indicates whether ju is 
an interval measure or complementary measure and whether or not ab > 0 or ab < 0. 
Then we wish to find the coordinates (u,v) of the children using the geometry in the 
uv-plane. By Lemma 3.4 we know that such a yz exists and is unique. 

Here are some thoughts before we dive into the technicalities of the algorithm. Let 
Ip = (Up, Vp) be the wv-coordinates of the parent interval [—1, 1], and let q: = (ui, v1) 
and q2 = (ug, v2) be the vectors of the children. The parent measure is always a convex 
combination of the children measures. Hence, q, will always lie on the line segment 


connecting q, and gz. The development of the algorithm often uses the connections 
X [r,s] 
(dz 


So at times we relax our vocabulary in order to get a clear picture of the behavior. 


between probability measure py of the form du = dx and the interval |r, s] itself. 


We will exploit the symmetry of the curves and relationships with respect to the 
v-axis. This is to be expected because of the symmetry in the transformations. For a 
probability measure ju, consider q/(j1) where q/(j1) is the reflection of q() across the 
line u = 0. Let q and qo be the children vectors found from the vector q(j1) and let qi} 
and q}, be there reflections. Then the children vectors for q’ are qi, and qs. Note that 
qi now corresponds to the right child and gq‘ corresponds to the left child. Thus for 
any measure 1, we find the children vectors from either q(j) or q/(4). Furthermore, 
the symmetry of all of the boundary curves and regions implies that q(w) and q‘({1) 
lie in the same region. Therefore, the following algorithms safely assume u, > 0. 

We first consider the simple case of u, = 0. Then gq, lies in either Ag, To, or 
(0,1/3). If q = (0,1/3), then the corresponding measure pz is of the form du = 
1/2dxz. Thus q, = (—1/2,1/3) and q = (1/2,1/3) after readjusting to probability 
measures. If g, € Ag, then pu is the probability measure corresponding to [-b,b]. 
The left child corresponds to [—b, 0] and the right child corresponds to [0, b]. These 


measures correspond to boundary case a = 0 or b = 0, which for A is the curve gj. 
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Combining the fact that uj = —u2, q@ and q lie on g;, and that q, is on the line 
segment connecting g, and q2, we have that gq, and q2 are the intersections of the 
constant line v = v, and the boundary curve g;. We have a similar argument for 
dq € YT. The representative measure corresponds to [—1, —b] U[b, 1], so the left child 
corresponds to [—1, —b] and the right child corresponds to [b, 1]. Measures of these 
forms correspond to the boundary curve gz. Thus q; and q2 will lie on the curve 
g2. Therefore, q; and q2 are the intersections of the constant line v = v, and the 
boundary curve gp. 

Consider the case g, € A;. Then the representative measure corresponds to [a, b] 
with ab < 0. This implies that the entire representative measure is contained within 
only one of the children. Since u, > 0 we have that gq; = (0,0) and q2 = q. 

When q, € Ag, the representative measure corresponds to [a, b] with ab < 0. Then 
the left child measure corresponds to [a,0] and the right corresponds to [0,}]. Just 
as we had before, this implies that g; and qo lie on the boundary curve g;. Consider 
qi = (ui, V1) and qo = (ug, v2) in terms of a and b. Equations (3.23) and (3.24) give 


us the following equations: 


a+0 a a+0-at+0 @ 
rr - 3 3 i286) 
O+b  b 0+0-b4+8? 6 
U2 5) 5) U2 3 3 (3 9) 
The line connecting q; and q2 has slope 
V2—V, _ 2(b?—-a?)  2%at+b) 4a+bdb_ 4 (3.30) 


i= BG=a) «Ok Be 

Note that in this case, the slope only depends on the u-value of the parent interval. 

Therefore, g; and q2 are the intersections of the line through q, with slope FUp and 
the curve gj. 

Next we consider g, € Y. Then we have that the representative measure ju corre- 


sponds to |[—1, a] U[b, 1] with a gap of length 2, = (b—a). Replacing b with a+ @, gives 


(2 


us that 4 corresponds to [—1, a] Ula+é,, 1]. First we find the gap length by considering 
the curve w given by the wv-vector of the measure corresponding to |—1, t] U[t+Z, 1] as 
t changes. Using equations (3.25) and (3.26) we have that the parametric equations 


for w are given by 


i = (t+ 0)? 
2(2— (t+ £) +t) 


2-(t+é% +8 


es 52 = 4 Oe): 


and v(t, @) = 


(3.31) 


By solving equation (3.31) for t in terms of u and then substituting into the expression 


for v, we have 
’ L 1 2u? 
é,{-1,1))= ee ye ee 3.32 
w(u, £, [-1, 1]) gt ae (3.32) 


Generalizing this for the measure corresponding to [r, t] U[t + ¢, s], we have that 


1 
(3r3 + 772 — 12r?u + 377s + 12u?r — 3rs? — 2rel — rl? 


w(u,f, |r, s|) = ioe 


+ 12u2é + 12s7u — 12u2s — 3s? + 572+ sé? + £°). 


Specifically, g, must lie on the curve v = w(u, @, [—1, 1]) for some ¢ = @,. Then we find 
é, by solving vp = w(Up, 4, |—1, 1]) for 2, which is equivalent to finding the appropriate 
zero of 


fe 2 1 
tee (u2 —U,)+ 5) € — Que (3.33) 


using Newton’s method. Note that @, € [0,2], f(0) = —2u?, and f(2) = 2 — 2u,. 
Furthermore, u € [—1, 1] implies that f(0) < 0 and v < 1 implies that f(2) > 0. By 
the Intermediate Value Theorem, there exists an ¢, € [0,2] such that f(¢,) = 0. 

When gq, € Y, the gap is entirely contained in one child. Since u, > 0, the 
right child corresponds to the whole interval [0,1] and the left child corresponds to 
[—1, a] U[b, 0] with the same gap length as the representative measure. Thus gq. = 
(1/2, 1/3). For the left child the gap length is known so q; is the intersection of the 
curve w(u, &,,[—1,0]) and the line through the points q, and q. 

The final case to consider is g, € Tz. The left child corresponds to [—1, a] and the 


right child corresponds to {b, 1]. So we have that q and q both lie on the boundary 
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curve gj and we also have that they lie on the line passing through q,. In a previous 
case, we showed there is a very nice relationship between the slope of the line and 
Up, but we do not have any such relationship here. Therefore, an iterative method is 
implemented to find the children vectors. The idea is to start with a value for u2 and 
then determine the corresponding gap length. Once the gap length of the iterative 
method matches the actual gap we have found the actual value of wz. From ue, then 
V2 = go(uz) and q, is found to be the intersection of the line through q2 and q, and 
the curve qo. 

The intersections of the curve gz and the line through the vector q, and the vector 
(0, w(0, 2, |[—1, 1])) serve as good initial guesses for ui; and uz, which we denote a 
and a0), The k*” step of the iteration is denoted a” and a”), We desire to have 
the sum of the gap in [—1,0] and in [0,1] equal to @,. Thus the quality of the 
approximation at step k is determined by the discrepancy between the gap length 
of the approximation and @,. Let @) be the gap length on [—1,0] and OW the gap 
length on [0,1] based upon the approximation a”). Fortunately, we do not need to 
solve w for the gap lengths as there is a simpler relationship. Instead we write the 
gap lengths as a function of a” and a”, We have that ¢, = b —a, @) = —d, and 


he = b. Furthermore, 
a” =1/2(@-1), and a =1/2(6+1). 
Thus we write @) = — 24”) —land OW) = 2a") — 1. The discrepancy at step k is 
d® = ¢, — (-20 —14 2a) —1) = 2, + 24 — 20” +2. 


Then we adjust by aft! = a* + d/2. The iteration terminates when d is suffi- 
ciently small. 

The algorithm so far has only addressed the vectors in wv. Let Qx({—1, 1], p) be 
the quantities supplied for some measure p. At this point, we need a way of approx- 


imating Qx([—-1, 0], e) and Qx((0, 1], ) with measures of the form du = Cw(x)dz. 
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Notice the addition of the constant C’. The vector corresponding to p is 


M,([-1, 1], 2) vane 
Mo([-1, 1], e)’ Mo([-1, 1], p) 


From the algorithm above we find approximations for the wv-vectors of the children 


q(p) = (Up, Up) = ( 


corresponding to an unknown representative measure ju of the form du = w(x)dx with 


(1 = Xa.) 


X [a,b 
w(x) Pl or 0] 0 ee ea 


~ b-a 


d 


such that g(p) = q(u). Thus 
i a ; w(x)dx and vp» = / : 2w(x)\dx 
Let C = M([-1,1],p). Then Mo([—1, 1], p) = C = Mo([-1, U, u). Furthermore, 
My((-1,1], 9) = upMoll-1, 1,0) =¢ f mole)de = Mi((-1, 1],4) 
and finally 
My({-1,1],p) = %pMo((-11], 9) =C f 2%u(a)de = My({-1, 1], 0) 


Thus Qx({—1, 1], e) = Qx([-1, 1], «). We know C but we do not explicitly know the 
parameters a and b for the measure yu. It remains to determine Mo({(—1,0], w) and 
Mo((0, 1], #). By definition Mo([—1, 1], 4) = Mo({[—1, 0], #) + Mo((0, 1], u). Then we 
write Mo([—1, 0], u) = tC and Mo((0, 1], uw) = (1 —t)C for some parameter t € [0, 1]. 
Think back to the wv-plane with q(p) = qd» and children q and q2. The closer the 
vector q; is to gp, the larger the weight placed on the corresponding child. Let d be 


the Euclidean distance in the uv-plane. Then we define the parameter t¢ as follows 


d(q2, dp) 
d(q, q2) 


Therefore, the children distribution-dependent quantities are 


d(q2,4p = A 42 dp 
Mo({-1,0], 4) = Ges}C —— Mo([0,1], 4) = (1 ~ Ge283,0) 


) 
M,([-1, 0], 4) = u1Mo([-1, 0], 4) Mi (0, 1], 4) = u2Mi([0, 1], 1) 


M3({—1, 0], #4) = v1Mo([-1, 0], u) ~Ma([0, 1], 4) = veMi((0, 1], 1). 
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Figure 3.8 The overlay of I from the polynomial weighted scheme and the regions 
from the subdomain scheme. The blending scheme combines the results from the 
two schemes into one scheme. 


Blending 


For a measure p with quantities Qx({—1,1],) we use the methods above to find 
an approximation for the children using either the polynomial weight or subdomain 
measures. The downside to the polynomial weight scheme is that there are measures 
that give quantities not possible for the polynomial scheme. That is to say the region 
I in wv corresponding to polynomial weight measures does not cover all of 2. 

The polynomial weight scheme can be extended to the regions outside of [’, but 
large errors will occur. To illustrate the errors that arise outside of I’, consider the 
point cloud D that contains 713 points evenly distributed along the interval [—12, 12] 
with a gap from —7 to —2, and the y-values y = 2x + 1. Figure 3.9 clearly displays 
the error that occurs around a gap when the polynomial weight scheme is used for 
any measure. 

On the other hand, all measures can be represented by a subdomain measure. 
When a measure is close to but not exactly Lebesgue, the representative measure 
found using the subdomain scheme assumes there is a gap. This is not good, because 


our point cloud applications will not have perfectly Lebesgue measure even when the 
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Figure 3.9 ‘The first image depicts the average y-values of the 6 initial intervals 
from which the actual quantities are calculated directly from the point cloud D. 
The other images are levels 1, 2 and 3 of subdivision where the quantities are 
calculated according to the polynomial scheme described in Section 3.1. These plots 
illustrate the issue that occurs when strictly using the polynomial scheme. Notice 
the large errors that occur around the gap. 


Lebesgue measure is the most appropriate. Therefore, we choose to blend the two 
schemes. Let Qx([—1,0],7) and Qx([0,1],7) be the quantities calculated from the 
polynomial scheme and Qx({—1,0], ~) and Qx((0, 1], ) be the quantities from the 
subdivision scheme. 


We approximate the real Qx with Qx for some parameter t¢ 


Qx([-1, 0], p) = (1 = t)Qx({-1, 0],7) = #Qx([-1, 0], 1) (3.34) 


Qx([0, 1], p) = (1 — t)Qx([0, 1,7) + tQx([0, 1], 4). (3.35) 


The polynomial approximation is not considered whenever q(p) lies outside of I. 
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Thus, we let ¢ = 1 for q(p) outside of I. Now assume that q(p) € T. The closer 
the measure is to Lebesgue, the larger the polynomial approximation weight. In the 
uv-plane, a measure being close to Lebesgue means that the vector q(p) is close to 


(0,1/3). Consider the line L(u) passing through q(p) and (0, 1/3), 


wale, 5 a 
L(u) = um) Ut 


Then L(u) intersects the boundary of I in two locations r and s. Let d be the standard 
Euclidean distance function. Without loss of generality assume that d(q(p),r) < 


d(q(p), 5). Then we define t to be the ratio of the distances 


d(q(p); (0, 1/3)) 
d(r, (0, 1/3)) 


For t > 1, q(p) is outside of the [ and we set t = 1. Notice that for q(p) on the 


— 


boundary of I’, we have that the weight of the polynomial approximation is zero. 


When q(p) = (0, 1/3) the weight for the polynomial function is one. 


3.2. APPROXIMATION OF THE VALUE-DEPENDENT QUANTITIES 


Now that we have the approximations for distribution-dependent quantities Q x ({—1, 0]) 
and @x({0,1]), it remains to find Qy({—1,0]) and Qy((0,1]). The value-dependent 
quantities are not as local as the distribution-dependent quantities. Therefore we in- 
clude the quantities of the neighbor intervals [—3,—1] and [1,3] in the calculations. 
Then we write the children quantities as a linear combination of the parent and 


neighbor quantities. 


Qy([-1,0]) = co@y([-3, -I]) + a@y([-1, 1) + @Q@y ([L 3) (3.36) 


Qy ([0, 1]) = doQy([-3, —1]) + diQy([-1, 1)) + d2Qy ([1, 3) (3.37) 


for some coefficients co, C1, C2,do,d,, and dz. The standard subdivision schemes re- 


produce polynomials, so we do the same here for the y-values. Hence, we use the 
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test functions 1,2 and 2x? to calculate Qy({[—1,0]) and Qy({0,1]). Equations (3.36), 


(3.37), and the test functions give us the following systems of equations: 


0 = 1 3 
ldpx / ldpx / ldpx | ldpx Co 
0 Pi ri 3 
[2px = ie xdpx [rox i, xdpx C1 
0 ss) 1 3 
/ x’dpx i x’dpx / a’dpx | Pdvx C2 
a -3 -1 1 
and 
; i cr i, “4d i “1d 
1 
0 ldpx 3 Px = Px i Px do 
| rdpx |= i rdpx 4, rdpx | rdpx dy 
A a ch 3 
[ dex i) x’dpx / x’*dpx | Pdox dy 
=: -1 1 


We already have an approximation for all of the integrals in these systems. So the 


systems simplify down to 


Co 


Qx((-1,0) = | Qx(-3,-1) | @x((-1,1) | @x(tt,3) | | es 


C2 


and 
do 


Qx((0.1) = | Qx(t-3,-m) | @x(t-1,1) | @x(t1,3) |] a 

dy 
We solve for the coefficients cg, C1, C2,do,d,, and dy. At this point we run into the 
same issue as Section 2.1. The systems could be near singular, causing difficulty 
directly solving. Therefore, we solve the systems just as before, using the truncated 
SVD for some threshold. Once the coefficients are found, we plug the values into the 
equations (3.36) and (3.37), giving us Qy ([—1, 0]) and Qy((0, 1]). From Qy and Qx an 
approximation of the surface can be constructed. Our focus is on the calculation and 


preservation of the quantities, not how the quantities are used in the approximation. 
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Numerical Results for Distribution-Dependent Subdivision 


over Intervals 


In this section we present the results of the blending algorithm over the intervals. 
The theory suggests that the algorithm should perform extremely well for quadratic 
functions and regular distributions with gaps. The domain X is the interval [-12,12] 
and the initial partition consists of evenly spaced intervals. The intervals are bisected 
at each level of subdivision. We insert the points from the point cloud into the 
intervals and directly calculate the distribution-dependent quantities Mo, M,, M2, and 
the value-dependent quantity M,. No other part of the algorithm directly accesses the 
point cloud. Let /* denote the partition where the quantities are calculated directly 
from the point cloud. The subsequent levels of subdivision calculate the quantities of 
the children based on the parent’s and parent’s neighbors’ quantities according to the 
blending algorithm in Section 3.1. The figures in this section colored magenta are the 
actual values (not averages). The green intervals are the actual average quantities 
calculated over the initial partition. Then the blue intervals correspond to the first 
level of subdivision, red is the second level, and black is the third level of subdivision. 
At each level we plot (M,/Mo), (Mi/Mo), and (M2/Mo) values, which correspond to 
the average y,x, and x”. 

For the first example, J* is six evenly spaced intervals partitioning X, and we 
consider the point cloud D, which consists of 626 points uniformly distributed from 


-12 to 12 with gaps [—9.213, —4.0234] and [3.1, 5.23], with y-values satisfying 
y(x) = —3a? + Qe +1. 


Indeed, the results for this situation are great (Figures 3.10 - 3.14). Subdividing to 
level 10 gives a good idea of the limit function of the subdivision scheme (Figure 
3.14). 


Expanding the type of functions, our second test function has that J* is 12 evenly 
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Figure 3.10 Graphs of the initial averages for y, x, and x? calculated from the 
actual point cloud Dy. 
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Figure 3.11 Level 1 of subdivision: Plots of the averages for y, x, and x? calculated 
only from the quantities of the parent and parent’s neighbors. 
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Figure 3.12 Level 2 of subdivision: Plots of the averages for y, x, and x? calculated 
only from the quantities of the parent and parent’s neighbors. 
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Figure 3.13 Level 3 of subdivision: Plots of the averages for y, x, and x? calculated 
only from the quantities of the parent and parent’s neighbors. 
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Figure 3.14 Comparision of the average y-value for level 10 of subdivision (blue) to 
the actual point cloud D, (magenta). 


spaced intervals over X, and the point cloud Dy, again consists of 626 points uniformly 


distributed from -12 to 12 with gaps [—9.213, —4.0234] and [3.1,5.23]. Now we let 


= sin x? 
a 18}. 


The algorithm performs well for point cloud Dg, although not as well as in D, (Figures 


the y-values satisfy 


3.15 - 3.19). One interval of interest is [8,10] € I*. Notice how the initial interval 
[8, 10] better approximates the local minimum as we refine strictly from the parent 


and neighbor values. 


3.3 APPROXIMATION OF Q OVER TRIANGLES 


Now we extend our ideas to d = 2 with triangles instead of intervals. We develop our 


algorithm for the simple case of 


Qx(R, w) = (Moo(R), Mio(R), Moi(R)) and Qy(R, uw) = My(R) (3.38) 


with 
Mo(R) =f tdu —— Myo(R) = [21 a (3.39) 
Mo(R)= fa du — My(R) = | y dp. (3.40) 


82 


Figure 3.15 Graphs of the initial averages for y, x, and x? of Dg. 
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Figure 3.16 Graphs of the averages for y, 2, and x? of Dz at level 1 of subdivision. 
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Figure 3.17 Graphs of the averages for y, 2, and x? of Dy» at level 2 of subdivision. 
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Figure 3.18 Graphs of the averages for y, 2, and x? of Dz at level 3 of subdivision. 
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Figure 3.19 Graphs of the average y for level 10 of subdivision (blue) and the 
actual point cloud D2 (magenta). 


Following the same ideas as the d = 1 situation, we make the assumption that the 
parent triangle is A* = ((0,1),(—1,0), (1,0)) and is divided into its left child A} = 
((0,0), (0,1), (—1,0)) and its right child A} = ((0,0), (1,0), (0,1)). This assumption 
is legitimate because there are affine transformations to and from a general right- 
isosceles triangle to the reference (standard) triangle A*. The transformations are 


defined here. Let (b,1,r) be a general right-isoceles triangle in X. Define T as 


—1 


i} 0 0 Le 
T= by ik —by ry — by 0 -1 1 


by ly — bo rg — be 1 -1 -l 


An element x* = (27, x3) of the reference triangle A* is transformed to x = (#1, 22) 


in the general triangle (b,1,r) by the affine transformation 


1 ib 
T a | ae 
5 a) 


Thus, elements from the general triangle is transformed to the reference triangle by 
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the inverse, 


1 1 
T*} 2, /= a 
L2 igs 
Define for a measure ps over A* 
Mio(A*, #) Mo (A*, 1) 
u(w) = ————— _ and v(p) :-= —————_. , 
4) Tag Be, uw) Me Teen B, 1) 


and q(u) = ((u(u), v(w)). We consider u and v as the average values of x, and x2 
respectively. Thus the vector q(j) will lie in the convex hull of A* for measures p. 
We develop a few methods based on the types of measures we wish to reproduce. The 


measures we consider are of the form dj: = w(x)dx with: 


e Weight Measures: 


W(X) = awa(x) + Buw(x) + YWwe(x) 


e Subdomain Measures: 


w(x)=Xe or w(x) =(1— xi) . 


Once the relationships in the uv-plane are established, then we introduce a constant 
C into the form for the measure. The weight measures have parameters (a, 3,7) and 
the subdomain measures have parameters based on the triangle t. The three parent 
quantities in Qx(A*) are used to determine the parameters. Once the parameters 
are found, then the children quantities are easily calculated. The last method on 
the triangle is to combine the two schemes into one blended scheme. The main idea 
on the triangle A* with Qx(A*) is to use the geometric behavior in the uv-plane 
to approximate Qx(Aj) and Qx(A3). To do this a one-to-one correspondence is 
established between the measure parameters and the appropriate region in the uv- 


plane. In this case it is possible to algebraically find the inverse. 


85 


Weight Measures 


Consider the probability measures pz over the reference triangle A* such that du = 
w(x)dx where w(x) = aw,(x) + Gw,(x) + ywe(x). We define the weight functions in 


such a way that 


is Wa(x)dxidry = i: wp(x)dx1dx2 = I. w-(x)drjdrg = 1. 


We also impose the following conditions on the weight functions 


w((0,1))=1 and w,((—1,0)) =w,((1,0)) =0 (3.41) 
w»((—1,0)) =1 and «,((0,1)) =ay((1,0)) =0 (3.42) 
w((1,0))=1 and w,((—-1,0)) =,((0,1)) = 0. (3.43) 


Therefore the weight functions are defined to be 


Wa(X) = Le (3.44) 
wp(x) = 1/2(—2%1 — z2 + 1) (3.45) 
w(x) = 1/2(4, — r2 +1). (3.46) 


Since yp is a probability measure, we have two conditions. The first condition is 
that w(x) is nonnegative for all x € A*. This implies more specifically that w is 


nonnegative at the vertices of A*. Then we have that 


w((0,1)) =a 20 (3.47) 
w((-1,0))=68 20 (3.48) 
w((1,0))=720 (3.49) 


The weight functions wz, Ww», and w, are all nonnegative so the inequalites above are 
sufficient to force w(x) > 0 for all x € A*. The second condition given by the fact 


that jz is a probability measure is that 


ie AWa(x) + Buy(x) + yw(x)dx=1 => (a +6+y7)=1. 
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Thus a = 3— 8—y making y < 3—§. Consider the transformation from the 
parameters of the weight measures (a,{,7) to the moments (Moo, Mio, Moi). We 
simplify this transformation to (3,7) to the moments (Mio, Moi) as Moo = 1 and a 


is a function of 6 and ¥. 


u= Mio = fay (3 — 8 —r)valx) + Bers(xe) + wel) dx 


= [6 —B-y)t1t_4+ PO an —Z2+1)+ Fle — 2+ 1)dx 


= 5 [23 8) + maa(6 - 36 - 37) + 21(8 + y)dx 


eaten (3.50) 


v= Mo = fea ((8 — B — y)walx) + Burl) + wel) dx 


7 f (8 — B 7) + (ay — 2 +1) + E20 — 22 + Ix 
i 
= = |. #3(6 — 37 - 38) + rr0a(y — 8) + 02(8 + y)dx 
2 2 12) (3.51) 


Therefore, 
(3,1) + (3(7- 8), (6-8 -7)) 
7 12 Y Tp Y) | 
This is a linear relation so the inverse is simply 3 = 3—6u — 6v and y = 3+ 6u — 6v. 


So we have 


(u,v) > (3 — 6u — 6, 3 + 6u — 6v). (3.52) 


The conditions on @ and ¥ are that 7,8 > 0 andy <3—. Then 6 > 0 implies that 
3 —6u —6v > 0, sov <1/2—u. The inequality 7 > 0 implies that 3 + 6u — 6u > 0, 
thus v < 1/2+u. Finally, the inequality y < 3— gives us that v > —1/4. Therefore, 
the image in the uv-plane is the triangle A,, with vertices (0, 1/2), (—1/4, 1/4), and 
(1/4,1/4). As to be expected, the weight measures are not enough to cover the range 


of all possible probability measures. 
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Figure 3.20 The graph is plotted in the wv plane, and the blue region is A,,. The 
triangle A,, indicates the vectors q for which a weighted measure corresponds. 


Here we outline the algorithm and introduce the constant C into the representative 


measure 7. Assume that p is a measure over A* such that q(p) € Ay. Then 


Mo(A", p) Moi (4%, p) 
a(p : Caron Tea) 


From q(p) we have 6 = 3 — 6u — 6v, y = 3 + 6u — 6v, and a = 12v — 3 from the 


inverse equations (3.52), such that 


=) xriw(x)dx and v(p) = Mal" 2) 


Let C = Moo(A%*, p). Then for the measure 7 with dy = Cw(x)dx, we have 


Moo(A", p) =C= Moo(d"*, n) (3.53) 
Myo(A*, p) = u(p)Moo(A*,p) = Cf ayw(x)dx = Muo(A*n) (3.54) 
Mo1(A*, p) = 0(p)Moo(A*, p) = C a _tao(x)dx = Moi(A", 7). (3.55) 


Therefore, from the unknown measure p quantities we construct the known represen- 
tative measure 7 such that Qx(A*, p) = Qx(A*,7). Next we directly calculate the 


quantities of the children over 7 to approximate the children quantities of p. So 


Qx(Aj, p) oo Qx (At, 7) and Qx(A3, p) Qx(A3, 7). 
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Figure 3.21 Illustration Figure 3.22 Illustration Figure 3.23 Illustration 
of a triangle measure of of a triangle measure of of a triangle measure of 
Type 1. Type 2. Type 3. 


Subdomain Measures 


In this section we develop an algorithm for representing an unknown measure p with 
a known measure y of the form du = C'y;dx or du = C(1 — yz)dx for some constant 
C’ and some triangle t. The first form is referred to as triangle measures and the other 
are complementary measures. Now we define the set D which specifies the types of 
triangles t we consider. Define three vectors that lie on the boundary of A* in the 


x-plane as follows: 
e b=(b,1—b), whereO0<b<1 
ec=(c,l+c), where —-1<c<0 
e d=(d,0), where —-1<d<l 


Let D denote the set of triangles of the form (b, d, (1,0)), (c,d, (—1, 0)), or (b, c, (0, 1)) 
for any b,c, and d in the appropriate ranges. Note the symmetry of the triangles in 
D with respect to the line x, = 0. The triangle (b,d, (1,0)) with d > 0 is the same 
triangle as (c,d, (—1,0)) with d < 0 reflected across the line x; = 0. In the same way 
(b,d, (1,0)) with d < 0 and (c,d, (—1,0)) with d > 0 are symmetric. Therefore, we 


have D separated into three types of triangles, disregarding the symmetric triangles 


e Type 1 (Fig 3.21): (b,d, (1,0)) with d>0 
Triangles of Type 1 are defined by the parameters b,d and C’, and are charac- 


terized as triangles contained completely in one child. 
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e Type 2 (Fig 3.22): (c,d, (—1,0)) with d>0 


Triangles of Type 2 are defined by the parameters c,d and C. 


e Type 3 (Fig 3.23): (b,c, (0, 1)) 


Triangles of Type 3 are defined by the parameters b,c and C’. 


Type 2 and 3 triangles intersect both children. 


Triangle Measures 


We begin by examining the transformations to the uv-plane for triangle measures p of 
the form du = C'y;dx over the different types of triangles t. This section determines 
the equations for the invertible relationship between the uv-plane and the parameters 
b,c and d, and defines the region of the wv-plane for which this relationship holds. 
We begin by noting that the u coordinate corresponds to the average x; value and 
the v coordinate corresponds to the average x2 value. So the image in the wv-plane 


will be contained within the convex hull of (-1,0), (1,0) and (0,1). 


Type 1 Triangle Measures _ For Type | triangle measures, the triangle t is bounded 


by the lines 


—2r,+d (3.56) 
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Thus the moments are 


Moo(A*, 2) = ee | Cie | ee ae i Canis (3.57) 
=f 1(b-1) 

My(A*, 2) = ie eae | Cpie= i a ne ee Cardade, (3.58) 
= = (d- 1)(b—1)(1 454d) 

Moi(A*, ) = ie ee | Cie = ae os Crodeydt, (3.59) 
7 <a —1)(b—1)". 


Therefore the transformation (b,d) — (u,v) is defined as follows 


My(A*,w) _1+b+d 


Mae) 


(3.60) 


and 
Moi(A*,) _ 1-6 
Moo(A*, |) 3 


The transformation is linear, so the inverse (u,v) — (b,d) is easily determined to be 


v(b, d) = 


(3.61) 


b(u,v) =1—3v and d(u,v) = —2+3v 4 3u. (3.62) 


Let R, denote the image of the square 0 <b <1 and0<d< 1 inthe uwv-plane. The 


Table 3.3. Type 1 Triangle Measures from the boundaries in the bd-plane to 
uv-plane 


Boundary in bd-plane | Transformation Image in wv-plane 
C0 w=1/311+d) v=1/3 eels 
ot u=1/3(2+d) v=0 v0 
a=0 w=1/3(1+6) v=1/31—)) | v=2/3-u 
ail u=1/3(2+6) v=1/31-—b) |v=1-u 


transformation (b,d) — (u,v) is a simple linear transformation, so the transformation 
of the boundaries of the square in the bd-plane are sufficient to describe R;. From the 
boundary calculations found in Table 3.3, we have that R, is the quadrilateral with 
vertices (1,0), (2/3, 0), (1/3, 1/3) and (2/3, 1/3) in the wv-plane. Therefore, measures 


corresponding to a vector gq € R,; can be modeled by a triangle measure of Type 1. 
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Type 2 Triangle Measures The triangle t for Type 2 triangle measures is bounded 


by 


v2 > 0. 
Thus the moments for Type 2 triangle measures are 


ct+1 Lea(x2) 
Moo(A*, 12) =f. du = | Céx = f ie Cdaidry 


C 
c+1 Lea(x2) 
M(A*, ws) ae. xydp = [ondx = | - Ca dzr,dx2 
= (c + 1)(1+d)(d+c-—1) 


ctl Lea(x2) 
Moi (A%, 1) =f. ndu= [ Cardx = | es Carodx, dx 


= = (c+ 1)?(1 +d) 


Thus we define the transformation (c,d) > (u,v) by 


-1 1 
oe and Cos 


aed) = 3 


The inverse relation (u,v) — (c,d) is defined by 


c(u,v) =3v—1 and d(u,v) =2—3v+4 3u. 


(3.63) 


(3.64) 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


(3.69) 


Let Rz denote the image of the square —1 < c < 2and0 <d< 1linthe wv-plane. The 


transformation (c,d) — (u,v) is a simple linear transformation, so the transformation 


of the boundaries of the square in the cd-plane are sufficient to describe Ry. From 


the boundary calculations found in Table 3.4, we have that R»2 is the quadrilateral 


with vertices (—2/3, 0), (—1/3,0), (0, 1/3) and (—1/3, 1/3) in the uv-plane. Therefore, 


measures corresponding to a vector qg € Ry can be modeled by a triangle measure of 


Type 2. 
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Table 3.4 Type 2 Triangle Measures from the boundaries in cd-plane to the 
uv-plane 


cd-plane | Transformation Image in wv-plane 
c=0 |u=1/3(-14+d) v=1/3 v= 1/3 

c=-1 |u=1/3(-2+d) v=0 v=0 

d=0 |u=1/3(e-1) v=1/38(14+c) | v=2/3+4u 
d=1) |w=¢/3 v=1/3(14+c) | v=u41/3 


Type 3 Triangle Measures Lastly, we consider Type 3 triangle measures. In this 


case, the triangle t is bounded by 


=l|- XY 
c+b)(a,—b 
— Lye(21) = ( Pes ) 1—- b. 
Thus the moments for Type 3 triangle measures are 
1+21 1-2 
Moo(A*, 12) =f : __ Cieraers + [ [. Odrade, (3.71) 
Loc £1) 
= —Cbc 
1+21 b 1-2 
Myo(A", 1) =e A Caydargda, + i | Caydrgdry (3.72) 
Lge (#1) 0 JLyc(x1) 
=“ be(b +c) 
1+21 1-2 
Mo:(A*, 12) a [ __ Cateraders + [ he Carndrdry (3.73) 
Loe (#1) Lpc(#1) 
=“ bel6 —c-—3) 
The transformation (b,c) — (u,v) is defined by 
b —b 3 
(D.C). = —- and v(b,c) = — (3.74) 
Then we have the inverse relation (u,v) — (c,d) which is 
b(u,v) =3/2.1—v+u) and c(u,v) =3/2(-l+v+4u). (3.75) 


From the transformation of the boundaries in the bc-plane, we describe the im- 
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Table 3.5 Type 3 Triangle Measures from bc-plane boundaries to the uv-plane 


bc-plane | Transformation Image in wv-plane 
b=0° | w=<c/3 v=1+c/3 v=l1+u 

b=1 |u=1/38(e+1) v=1/38(24+c) | v=1/3+4u 
c=0 |x=0/3 v=1-0)/3 v=l1-u 
ec=-1 | u=1/3(6-1) v=1/3(2-—b) | v=1/3-u 


Figure 3.24 The uv-plane colored according to the types. Type 1 is blue and the 
symmetric case is light blue. Type 2 is purple and the symmetric case is light 
purple. Type 3 is yellow. 


age of the square 0 < b < 1 and —1 < c < O in the we-plane. Table 3.5 con- 
tains the boundary calculations. Hence, the image is the square R3 with vertices 


(—1/3, 2/3), (0,1), (0,1/3) and (1/3, 2/3). 


Table 3.6 The Parametric Equations and Inverses for the Triangle Measures 


Type | Transformation Inverse 
1+b6+4d T—b 

1 ue aes b=1-3v d = —24+ 3u+ 3u 

—1 1 

2 w= tet v= 6] a0 = d = 2—3v+3u 
b+c —b+c+3 3(1-—v+u) 3(-—l+u+u) 

3 U= nn b= CS SS 
3 3 2 2 


For each type of triangle measure, equations have been presented for the calcu- 
lation of the parameters b,c, and d given q = (u,v) of the measure over the parent. 
These parameters determine the triangle t. In the case of intervals, the constant C’ 
was simply the zeroth moment of the parent. In the case of triangles, it is not that 


simple. 
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Lemma 3.6. Let p be any measure over A* with q(p) in Ri UR2U R3 (or in the 
reflection of Ri UR2U Rs across the v-axris). Then there is a unique 4s of the form 
du = Cyx:dx with t € D such that Qx(A*, uw) = Qx(A", p). 


Proof. First identify in which region q(p) = (u(p),v(p)) lies. If q(p) € R; for i € 
{1,2,3}, then the representative measure is of Type 7. 


For i = 1, using the inverse relationship Eq (3.62) let 


b=1—3v(p) and d= —2+4 3v(p) + 3u(p). (3.76) 
Then let 
= 2Moo(A*, p) 
CL = (@-—1) (3.77) 


Define t,q to be the triangle with vertices (b, 1 — 6), (d,0), and (1,0) and py to be the 


measure such that du = Cx1,,dx. Equations (3.57), (3.58), and (3.59) then give us 


that 
ar Gare cs — — Moo(A*, p) 
Moo(A"*, 4) = 3 4 =)b=1)= d-6-)"” 1)\(o=1) 
= Moo(A", p) 
Mio(A*, 1) = Sd — 16-1) tb+d) = Moot") + b +d) 
2 Moot ) 4 +1 —30(p) —2-+30(p) + 3ulp)) 
= Moo(A", p)u(e) = Mio(A*, p) 
Mox(A*, 1) = =a —1)(0 — 1)? = OD) 1) 
= Moo(A", p)v(e) = Moi(A", p) 
Thus Qx(A*, 4) = Qx (4%, p). 
For 2 = 2, the approach is similar. Let 
c= 3u(p)—1 and d=2-— 3v(p) + 3u(p), (3.78) 
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as in equation (3.69). Then let 


= 2Moo(A*, p) 
aS (d+1)(c +1) om 


Define t.qg to be the triangle with vertices (c, 1+ c),(d,0), and (—1,0) and yz to be 


the measure such that du = C2xz,,dx. Equations (3.64), (3.65), and (3.66) then give 


us that 
Moo(A*, #) = S2(d + 1)(e-+ 1) = Moo(A*,p) 
Mi(A*, 2) = = (e+ 1)(d+e—-1) = Moot Ph Le-1) 
= Moot") (2 — 80(p) + 3u(p) + 30(p) — 1-1) 
= Moo(A*, p)u(p) = Mio(A*, p) 
Mu (At, 2) = (d+ I(e+1?2 = Moot?) (c+1) 


= Moo(A*, p)v(e) = Moi(O", p). 


Thus Qx(d"*, LL) = Qx (A*, p): 
For i = 3, let 


b= 3/2(1—v(p)+u(p)) and c=3/2(-1+4+ v(p)+u(p)), (3.80) 


as in equation (3.75). Then let 


— Moo(A*, 
C3 = a (3.81) 


Define t,- to be the triangle with vertices (b, 1 — b),(c,1+ cc), and (0,1) and yz to be 


the measure such that du = C3x,,,dx. Equations (3.71), (3.72), and (3.73) then give 
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us that 


Moo(A*, w) = —C3be = Moo(A", p) 


Myo(A*, n) = pels te) = Moot Py +e) 
= Moot") (3/2(1 — w(o) + u(p)) + 3/2(—1 + v(0) + u(o)) 
= Moo(A*, p)u(p) = Mio(A*, p) 
Moi (A*, 1) = be(—b be+3)= tt He+3) 
_ — Mel") (3/210 — u(p) + u(p)) + 3/2(-1 + v(p) + u(p)) +3) 


= Moo(A", p)v(e) = Moi(A*, p)- 


Thus Qx(A*, 1) = Qx(A*, p). 

It remains to find yz for q(p) in the reflection of Ri U R2U R3 across the v-axis. 
Suppose q(p) is in the reflection of some R; for i € {1, 2,3}. Then we replace q(p) = 
(u,v) with q/(p) = (—u,v). From the above work, we calculate the measure p’ of 
the form du’ = Cy,dx with q(u') = q'(p). Let ww be the measure such that du = 
Cy:d(—2x1)dx. Define t’ as the triangle t reflected across the x2-axis. Then du = 
Cxwdx. Thus Qx(A*,u) = Qx(A%p): 


Complementary Measures 


We turn our attention now to the complementary measures p° of the form du’ = 
C(1 — x,)dx (Figures 3.25, 3.26, 3.27). We start by examining the transformations 
to the wv-plane for each type. The complementary measures do not have the simple 
linear relationship with the wv-plane exhibited by the triangle measures. Thus we 
first show for each type that there exists a unique transformation from the wv-plane 
to the parameters b,c, and d. Later we show the calculation of this transformation. 


For the computation of the moments for complementary measures, we us the 
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So 


Figure 3.25 Illustration Figure 3.26 Illustration Figure 3.27 Illustration 
of a complementary of a complementary of a complementary 
measure of Type 1. measure of Type 2. measure of Type 3. 


Lebesgue measure over the entire triangle A*. Note that 


1 1-22 

ldx = - 7 ldx,dxy = 1 (3.82) 
A* 0 x2—-1 
1 1—2x2 

| 10x = | ; x dx ,dxr_ = 0 (3.83) 
A* 0 Jaxo—1 
if 1-29 

| todx = | / rodxidry = 1/3 (3.84) 
A* 0 Jao—1 


and Lebesgue measure over the entire triangle A* corresponds to the wv-vector 


(0,1/3). In future considerations we need not compute the transformations for 


(0, 1/3). 


Type 1 Complementary Measures The moments in terms of b,d, and C for the 


complementary measures of Type 1 are 
Moo(A*, u°) = i 1dye ={.. Cax — | Cax (3.85) 
= C —C/2(d— 1)(b—-1) 
Myo(A*, nw) = ie due = I. Ca,dx — [cedx (3.86) 
= —C/6(d—1)(b—1)(1+b+d) 
Moi(A*, 2") = I. rod = I. Caodx — [ Coodx (3.87) 
= C/3+C/6(d—1)(b-1)° 


Let M® be the transformation (b,d) + (u,v) given by 


— My(A*, nu) — (d-1)(b—1)(1+b+4) 
ue) = Moo(A*, 5) -6 +. 3(d —1)(6— 1) co) 
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and 


— MolAt2) _ (d@-1)6—1)2 42 
M8) = ar (Ae, 2) 6 —3(d—1)(b—1) eee) 


ME is continuous for all 6 and d in [0,1]. Thus examination of the transformation of 
the boundaries is sufficient to determine the region of the image. M. : of the boundaries 


of the square 0 < b < 1 and 0 < d < 1 from the bd-plane to the uv-plane are given 


2—9u? 
in Table 3.7. Hence, the image is bounded by the curves v = see le i= 1/3. 
3(2 + 3u) 
Table 3.7 Type 1 Complementary Measures from bd boundaries to wv 
bd-plane | Transformation Image in wv-plane 
b= 0" | wat sid). -e S173 w= 1/3 
b= 1 | | |= 0 oS 193 (0, 1/3) 
1+2b-06? 2— 9u? 
d=0 = 1/3(b-1 = = > TN 
UMM, De sea i 2 sorean) 
d=” || 4620 v=1/3 (0, 1/3) 


Note that the uv-vector (0,1/3) corresponds to Lebesgue measure over the entire 


triangle A*. Define 5; to be the set {(b,d) : b,d € [0,1)}, and let 7, be the region 
2—9u? 


3(2 + 3u) and v = 1/3 with a hole at 


in the wv-plane bounded by the curves v = 


(0, 1/3). 
Lemma 3.7. The map Me : 5S; 3 T, ts invertible. 


Proof. Note that M®(S,) C T, from the boundary analysis above. The Jacobian of 
Me is 
(d — 1)(b— 1) [(d — 1)?b? — 2d?b — 2db + (d+ 1)?| 

9(2 — (d— 1)(6— 1))° 


We wish to show that J # 0 for all (b,d) € S,. The factors (d — 1) and (b — 1) 


J= 


are disregarded as b # 1 and d # 1. Dissecting and regrouping the terms of the 


numerator, we have that 


(d — 1)*b? + (1 — d*b) + (d? — d?b) + (2d — 2db). 
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Each term is greater than or equal to zero as b,d € [0,1). Furthermore, the term 


(1 — d’b) is strictly greater than zero because b,d € [0,1). Thus the Jacobian of M? 


is nonzero on 5S). Hence we have that Mt is invertible 


Type 2 Complementary Measures We start by defining the map ME. Using 


the relationship between the complementary measures and the triangle measures, we 


have that 
Moo(A*, po) = C — C/2(d + 1)(c +1) (3.90) 
Myo(A*, pw?) = —C/6(d + 1)(c + 1)(c +d -1) (3.91) 
Mor(A*, pw) = C/3 — C/6(d + 1)(c + 1)?. (3.92) 


Therefore, the map Ms is given by 


_ (dt (c+ (e+d-1) _ (dt 1j(e+1)?-2 
Uae iichla6 | Gan iea =e ee) 


M§& is continuous for (c,d) 4 (0,1). Thus examination of the transformation of the 
boundaries is sufficient to determine the region of the image. Ms of the boundaries 
of the square —1 <c <0 and 0 <d <1 from the cd-plane to the wv-plane are given 


in Table 3.8. We have the additional boundary v = 1 — u, since u corresponds to the 


Table 3.8. Type 2 Complementary Measures from cd boundaries to uv 


cd-plane | Transformation Image in wv-plane 
6=0 | “= T3140) w]173 ves 
c=-1 |u=0 oH 178 v= (0) 173) 
—14+2c+¢e¢ Qu? — 2 
d=0 |u=1/3(1+ 22S aS 
Be ee)! Baa cay || aga) 
d=1 |u=1/314+c) v=1/38(2+c) v=ut4+1/3 


average x, and v corresponds to the average xr. This boundary also appears as a 


limit when both c + 0 and d — 1 and is discussed in Section 3.3. Hence, the image 
9u? — 2 


TG TSE ice ea and v= 1—-—u. Let 


is bounded by the curves v = 
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Sy = {(c,d) € (—1,0] x [0,1]\(0,1)}, and let 74 be the region in the uv-plane with 
Qu? — 2 


2 sag P= Vay Sutisv<i—4, and a hole at (0, 1/3). 


Lemma 3.8. The map Ms : Sy > T is invertible. 


Proof. The Jacobian of M$ is 


(d+1)(c+ 1) [(d+1)?c? + 2d?c — 2de + (d — 1)?] 


= | 
a= 9(2 — (d+ 1)(c+1))3 


We show that the Jacobian is nonzero for all (c,d) € Sj. The factors (d+ 1) and 
(c+ 1) are disregarded as c £ —1 and d #4 —1. Next we examine the rest of the 


numerator, which is 


(d+ 1)?c* + 2d(d — 1)c¢+ (d—1)?. 


Each term is greater than or equal to 0 since c € (—1,0] and d € [0,1]. If c £0, then 
the term (d+1)?c? > 0. Ifd 4 1, then the term (d—1)? > 0. Note that (0,1) ¢ S2, so 


the numerator of J is strictly greater than zero. Thus, the Jacobian of Ms is nonzero 


on S». So we have that Ms is invertible. 


Now we repeat the process above and define the map M$. We have that 


Moo(A*, p°) = C + Che (3.94) 
Myo(A*, p") = C/3be(c + d) (3.95) 
Moi (A*, uw") = C/3 — C/3bc(b — c — 3) (3.96) 


Then the map M$ is given by 


_ bc(c +b) oa . Seb e= 3) 1 


E30 4b) SSC) io) 


The map M$ is continuous for all (b,c) € [0,1] x [—1,0]\(1,—1). Table 3.9 contains 
the boundary calculations. We have the additional boundary v = 0, since u corre- 
sponds to the average x7, and v corresponds to the average x2. This boundary also 


appears as a limit when both b > 1 and c + —1 and is discussed in Section 3.3. 
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Table 3.9 Type 3 Complementary Measures from bc boundaries to uv 


bc-plane | Transformation Image in wv-plane 
b=0 |u=0 w= Ts (0,173) 

b= I> | |S e73 v=1/3(14+c) | v=1/3+u 
C0) |= O O= 1/3 (0, 1/3) 

c=-1 |u=)/3 v=1/3(1-—6) | v=1/3-u 


Hence, the image is bounded by the curves v = 1/3 —u, v = 1/3+u, and v = 0. Let 
S3 = {(b,c) € (0, 1] x [-1,0)\(1, —1)}, and let 73 be the region in the wv-plane with 
v<1/3-—u,v <1/3+ 4, v > 0, and a hole at (0, 1/3). 


Lemma 3.9. The map M$ : S3 — T3 is invertible. 


Proof. The Jacobian of M$ is 


2bc(b?c? + 3bc + b — c) 
J a: 
9(1+ bc)? 


Next we show that the Jacobian is nonzero for all (b,c) € S3. The factor 2bc is 
disregarded as b 0 andc 4 0. Then the numerator is rewritten with the substitution 


c = —c so that now c € (0, 1] instead of in [—1,0), and denote it by N 
N =06'c? —3bc+b+¢, (3.98) 


and N 4 0 if and only if J £0. Then 6 and c nonzero allows us to divide by bc. So 


The terms bc,1/c, and 1/b are all nonnegative so we apply the inequality of the 
arithmetic and the geometric mean. We use the strictly greater than form since we 


do not have that all the terms are equal. This gives us that 


~ > 3ifbe-1/c-1/6-3=3-3=0. 


Thus the Jacobian of M$ is nonzero on $3, and we have that M$ is invertible. 
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Figure 3.28 The uv-plane colored according to the triangle and complementary 
measure types. R; are triangle measures and T; are complementary measures for 
i= 1,2,3. 


Limiting Cases on the Boundary 


Lemmas 3.7, 3.8, 3.9 show that the relations are invertible for the vector gq in A* 
except for (0,1/3) and some of the boundaries. The uv-vector (0,1/3) is associated 
with the Lebesgue measure over A*, so no further analysis is needed for this vector. 
There are three sections of the boundary not included in the regions R;, T;, and their 


symmetric counterparts: 
1. v=1-—-u for u € (1/8, 2/3) 
2. v=1+u for u € (—2/3, 1/3) 
3. v =0 for u € [—1/3, 1/3]. 


These boundary sections correspond to limiting cases which appear when the area of 
the subdomain defining the complementary measure (A* t) is arbitrarily small. In 
these situations, the measure is concentrated along an edge. The first arises as a lim- 


iting case for Type 2 complementary measures. The second arises for the symmetric 
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case for Type 2 complementary measures. Lastly, the third arises for a limiting case 
for Type 3 complementary measures. 

The Type 2 complementary measures limiting case occurs when c = 0 and d= 1. 
Notice the ‘and’ statement here. This situation is similar to the limiting case we 
managed with the intervals. Depending on the way (c,d) approaches (0,1), we will 


have different uv-coordinates. Start by considering c = 0 and d approaching 1. Then 
u=1/3(1+d) 52/3 and v=1/3. 
On the other hand when d = 1 and c approaches 0, then 
u=1/3(1+c) 31/3 and v=1/3(24+c) > 2/3. 


Let 7 > 0 and € > 0, consider c = € and d= 1-—vTe. Ase > 0 then c > 0 and 


d— 1. We rewrite u and v interms of 7 and «. Then € — 0 implies that u approaches 
2(7 — 1) 
3(7 — 2) 
in (1/3, 2/3) since 


and and v approaches . Note that v = 1—w and wu can be any value 


eet 
3(7 — 2) 
lim u =1/3 and lim u= 2/3. 

Thus this limiting case covers v = 1—u for u € [1/3, 2/3] in the wv-plane. In this case, 
the measure is concentrated on the edge v = 1— u. Therefore, we consider the left 
child A} to be empty and the quantities of the right child Qx(A3) = Qx(A*). The 
corresponding symmetric case covers v = 1+ u for u € [—2/3,1/3] and we consider 
the quantities of the left child Qx (Aj) = Qx(A*) and the right child A$ to be empty. 

The limiting case in complementary measures of Type 3 occur when b = 1 and 


c = —1. First consider c = —1 and b approaching 1. Then 
u=6/331/3 and v=1/3(1-)) 50. 

On the other hand when b = 1 and c approaches 0, then 
u=c/3>4-1/3 and v=1/3(1+c) 0. 
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Now for 7,¢€ > 0, consider b = 1—e€ and c= —1+Te. As € approaches 0, u approaches 
T-1 


—_§— andv-— 0. Notice that u can be any value in (—1/3,1/3) since 
—3(7 +1) 


lmu=1/3 and lim u=-—-1/3. 
T—>0 T—+00 


Hence this limiting case covers v = 0 for wu € [—1/3,1/3] on the uv-plane. In this 
case, the measure is concentrated on the edge v = 0. If u > 0, then we consider 
Aj to be empty and Qx(A3) = Qx(A*). If u < 0, then we consider A} to be 
empty and Qx(Aj) = Qx(A*). If u = 0, then we let Qx(Aj) = $Qx(A*) and 
Qx(A3) = 3Qx(A*). 


Computation of Parameters for Complementary Measures 


Given the quantities Qx(A*,) for some measure p we wish to find a representa- 
tive measure such that Qx(A*,u) = Qx(A*,p). In Lemma 3.6 we have given a 
method for determining the parameters for triangle measures (q(p) € Ri UR2U R3 
or symmetric). Now we develop the method for complementary measures where 
q(p) € T, UT2UT3 or symmetric, in which case the representative measure is of the 
form p° with du’ = C(1 — x,)dx. The idea is to write the parameters b, c,d in terms 
of C and then use them to solve for C’. This requires solving a quadratic equation. 
Once we know the parameters and C, then the representative measure is fully known 
and the children’s distribution-dependent quantities can easily be calculated. We will 
need to look at each type individually. 

For Type 1 complementary measures we have the following system of equations 


that we need to solve for (b,d,C) given the moments of t°. 


C 


a (4 —1)(-—1) =C— Moo (A*, 0°) (3.99) 
= (d- 1)(b—1)(1 +b + d) = —Myo (A*, 0°) (3.100) 
<a —1)(b= 1)? =C/3 — Mu (A*,p") (3.101) 
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Dividing (3.101) by (3.99), we have that 


_ C/3 = Mor (A*, p®) 


—(6-1)/3 = 3.102 
—C + 3M (A*, p°) 
Sib + 1 
CS Moo(A*, p*) 
Dividing (3.100) by (3.99) and substituting (3.104), we have that 
—M,,(A*. of 
1/3(1 +b +4) 1o(A*, #”) (3.103) 


~ C= Moo(A*, 2°) 
7 —3Mi(A*, p*) + C — 3M (A*, p°) _ 


2 
Cs Moo(A*, p) 


>d 


Now we plug 6 and d back into (3.99) and solve for C’, which requires solving the 


quadratic function 
3C? (Mio — Mor + Moo) + 3C (—3MioMor — 3MG; + 3MooMor — 2Mgy) + 2Mgy = 0 


We need to be sure that Myo + Moo — Moy is not zero, which occurs when v = u + 1. 
Fortunately, the complementary Type 1 is nowhere near satisfying that equation. We 
reject the solution for C which implies that 6 or d is outside of their ranges. 

Using the same strategy as in Type 1, we have for complementary measures of 


Type 2 that 


_C/3= Moi(A*, p°) 
C= Moo(A*, p*) 

CPs 3Moi(A*, p*) 
C= Moo(A*, p°) 


1/3(c + 1) (3.104) 


>C —1 


Dividing (3.100) by (3.99) and substituting (3.104), we have that 


— —Mio(A*, p°) 

C= Moo(A*, p*) 

ater ae —3Myo(A*, p°) — C + 3Moi(A*, p*) 
C — Moo(A*, p*) 


1/8(etd =1) (3.105) 


2 
Now we plug c and d and solve for C’, which requires solving the quadratic function 
3C? (—Mio — Moi + Moo) + 3C (3MioMor — 3Mo; + 3MooMor — 2Moo) + 2Mop = 0 
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We need to be sure that —Mj9 + Moo — Mo is not zero, which occurs when v = 1—u. 
This only occurs in the limiting cases which we have already discussed in Section 3.3. 
We reject the solution for C’ which implies that c or d is outside of their ranges. 


Again we use the same strategy as in Type 1. 


_ Of/8 — Mu(A*, 9) 


—1/3(b-—c-—3)= 3.106 
—C + 3Moi(A*, p°) 
>b-—c= + 3 
C= Moo(d*, p°) 
We also have that 
—~3My9(A*, p’) 
b+c= 3.107 
C— Moo(d*, p°) ( 
Adding the equations (3.106) and (3.107) we have that 


2C — 2Moo(A*, p*) 


Then 


=3Migl(A* pC =3Mei(A*p') 
sen ee 3.109 
3 IC ae 2Moo(A*, pe) / ( ) 


Now we plug in } and ¢ and solve for C, which requires solving the quadratic 


function 
1207Mo +C (-9Mz, + 9M5, — 18Mo0Mo — 3Moo) + 4Moq = 0 


We need to be sure that Mo; is not zero, which occurs when v = 0. The complemen- 
tary measures of Type 3 have a limiting case there which we previously discussed in 
Section 3.3. We reject the solution for C which implies that b or c is outside of their 


ranges. 


Blending 


Let p be a measure with quantities Qx({—1, 1], e). Then we use the methods above to 


find an approximation for the children using either the weight or subdomain measure 
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Figure 3.29 The uv-plane highlighting the overlap of the weighted measures and 
the subdomain measures. 


schemes. The downside to the weight scheme is that there are measures that give 


quantities not possible for the weighted scheme. That is to say triangle 


Aw = ((0,1/2), (—1/4, 1/4), (1/4, 1/4)) 


does not cover the whole triangle ((0,1), (—1, 0), (1,0)). On the other hand, all mea- 
sures can be represented by a subdomain measure. When a measure is close to but not 
exactly Lebesgue, the representative measure assumes there is a gap in the measure. 
Similar to the interval case, this is not good because our point cloud applications will 
not have perfectly Lebesgue measure even when the Lebesgue measure is the most ap- 
propriate. Therefore, we choose to blend two schemes. Let Qx (Aj, 7) and Qx (A$, 7) 
be the quantities calculated from the weighted scheme and Qx (Aj, 4) and Qx(A53, 1) 
be the quantities from the subdomain scheme. There is no weighted approximation 
whenever q(p) lies outside of A,,. In this case, we use only the subdomain measure 
results. 

Now assume that g(p) € A,,. Let (Xo, A1, Az) be the barycentric coordinates of q(p) 
with respect to the triangle A,,, and denote the minimum by \,, = min{Ap, Ai, A2}. 
When q(p) is on the boundary of A,,, the smallest barycentric coordinate would 
be zero (Am = 0). When q(p) = (0,1/3), the vector corresponding to the Lebesgue 


measure, we have that all of the barycentric coordinates are equal to 1/3, so Am = 1/3. 
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Therefore, we approximate the real Qy with Qx given by 


Qx(Aj) = 3AmQx(Aj,7) + (1 — 3Am)Qx(A;, 1) (3.110) 


ford. = 1.2. 


Approximation of the Value-Dependent Quantities 


Now that we have the distribution-dependent quantities Q x (Aj, dpx) and Qx(A5, dpx), 
it remains to find Qy(Aj,dpx) and Qy(A3,dpx). The value-dependent quantity is 
not as local as the distribution-dependent quantities. Therefore we include the quan- 
tities of the three neighbor triangles Ag, Az, and Ap in the calculations. Here 
we use a fairly straight forward method to calculate the value-dependent quantities 
of the children as a linear combination of value-dependent quantities of the parent 
and parent’s neighbors. Let Ag, Az, and Apr be the neighbors of A* of the same 
level. More specifically, Ag shares the edge (—1,0) and (1,0), A; shares the edge 
(—1,0) and (0,1), and Ag shares the edge (0,1) and (1,0). If desired, one could 
substitute another method for solving for Qy(Aj,dpx) and Qy(A3,dpx) as long as 
it depends only upon the distribution-dependent quantities of the parent (A*), the 
parent’s neighbors (Ag, A;, Ag), and the children (Aj, A3). 


Then 
Qy (Aj) = @Q@y(Ag) + aQy(Az) + eQ@y(Ar) + c3Qy(A") (3.111) 
Qy(A3) = doQy (Az) + diQy(Az) + d2Qy (Ar) + d3Qy(A*) (3.112) 


for some coefficients Co, C1, C2, €3, dp, dj, dz and dg. We also desire the approximation 
to preserve the polynomial y-values. Hence, we use the test functions 1,7, and x2 to 
calculate Qy (Aj) and Qy(A3). That leaves an additional degree of freedom, so we 
add the requirement that the weights for the left and right triangles are symmetric. 


The idea is that the weights should be able to balance. Equations (3.111), (3.112), 
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and the test functions give us the following system of equations: 


ldpx | ld i ld i ld | ld 
Vi — PX = Px oe Ox! fe nUPx co 
he x \dpx _ des x\dpx Je x \dpx [, dex a xyidpx C4 
fo eed | x2dpx | x2dpx i, xodpx J e2dpx C2 
1 B L R 
0 0 1 iD 0 es 
and 
i¢ ldpx | ldpx ldpx i, ldpx ‘i ldpx a 
2 Ap Ar AR A* 
J erdex i, xdpx | rdpx | rdpx | x\dpx dy 
AS = Ap AL AR A* 
| r2dpx i rodpx | rodpx | rodpx | rodpx dy 
AS Ae Re, Ap A* 
0 0 1 1 0 ds 


We already have an approximation for all of the integrals in these systems. So the 


systems simplify down to 


Qx(Aj) Qx(Azp) | Qx(Az) | Qx(Ar) | Qx(A*) C1 
0 ° 0 1 1 0 C2 
and 
do 
[axa |_| @x(ds) | x(42) |@x(Ax) [@x(ar) | | a 
0 0 i 1 0 dy 
d3 


We solve for the coefficients cg, C1, C2, €3, dp, dy, dy and d3. At this point we run into 
the same issue as Section 2.1. The systems could be near singular, causing difficulty 
directly solving. Therefore, we solve the systems just as before, using the truncated 
SVD for some threshold. Once the coefficients are found, we plug the values into the 


equations (3.111) and (3.112), giving us Qy (Aj) and Qy(A3$). From Qy and Qx an 
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approximation of the surface can be constructed. Our focus is on the calculation and 


preservation of the quantities, not how the quantities are used in the approximation. 


Numerical Results for Distribution-Dependent Subdivision 


over Triangles 


In this section we present the results of the blending algorithm over the triangles. 
The theory suggests that the algorithm should perform well for linear functions and 
regular distributions with gaps of the forms above. For the following tests, the domain 
X is a 13 x 13 square centered at the origin. The initial partition consists of two 
right isosceles triangles and then is subdivided according to newest vertex bisection. 
After four levels of uniform subdivision (32 triangles), we insert the points from 
the point cloud into the triangles and directly calculate the distribution-dependent 
quantities Moo, Mio, Moi, and the value-dependent quantity M,. No other part of 
the algorithm directly accesses the point cloud. Let 7* denote the triangulation 
where the quantities are calculated directly from the point cloud. The subsequent 
levels of subdivision calculate the quantities of the children based on the parent’s 
and parent’s neighbors’ quantities according to the blending algorithm in Section 3.3. 
The figures in this section colored red are plots of the average y-value (M,,/Moo) over 
the triangles. The green figures are the average x,-values (Mjo/Moo), and the blue 
figures are the average x2-values (Mo, /Moo). 

Our first test is on the point cloud D;, which consists of 653,629 points uniformly 
distributed over X with a gap between the lines x2 = —0.97x, + 0.83 and x2 = 


0.17, — 1.706. The y-values satisfy the equation 
y(%1, 22) = 2(a1 — 1)(a1 + 2) — 8x0(x2 + 0.5). 


Figure 3.30 displays the actual point cloud D, with the triangulation 7*. Figure 3.31 


is the average values of the triangles for the triangulation 7* which is the basis for the 
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Figure 3.30 The point cloud D, and the triangulation 7*. Note that measure over 
each can be approximated well by the triangle and complementary representative 
measures discussed in Section 3.3. 


Figure 3.31 The average values over the triangles in the triangulation 7* and for 
the point cloud Dj. 


subsequent subdivisions. The three levels of subdivision of 7* are shown in Figure 
3.32 in which the quantities are calculated from the previous level. As we progress 
through the different levels of subdivision the algorithm approximates the average y 
well and the average 7; and x2 extremely well. The gap becomes increasingly apparent 
with each level of subdivision. Continuing the subdivision to level 10, we get the idea 
of the limiting surface (Figure 3.33). When comparing the algorithm results Figure 
3.33 and the actual point cloud Figure 3.34 we see that the algorithm result is very 
close to the actual point cloud with the exception of small jagged triangles along the 
edge of the gap. 

The second test is on the point cloud D2, which consists of 376,451 points uni- 


formly distributed over the circle centered at the origin with radius 5. The y-values 
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Figure 3.32. The rows depict the subdivision levels 1 through 3 of 7*. The red 
plots are the average y-value. The green plots are the average x;-values and the 
blue plots are the average x2-values. 


satisfy the linear equation 


y(@1,%2) =14+ 21 — 2. 


Figure 3.35 displays the actual point cloud D2 with the triangulation 7*. Figure 
3.36 is the average values of the triangles for the triangulation 7* which is the basis 
for the subsequent subdivisions. The three levels of subdivision of 7* are shown in 
Figure 3.37 in which the quantities are calculated from the previous level. As we 
progress through the different levels of subdivision the algorithm approximates the 


average y well and the average x; and x2 extremely well. Continuing the subdivision 


113 


Figure 3.34 The point cloud D, colored red according to the average y-value, 
colored green according to the average x,-value, and then colored blue according to 
the average £2-value. 


Figure 3.35 The point cloud D2 and the triangulation 7*. The better these initial 
measures over the triangles are represented by the triangle and complementary 
measures discussed in Section 3.3, the better the limiting surface of the algorithm. 
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Figure 3.36 The average values over the triangles in the triangulation 7* and for 
the point cloud D3. 


Figure 3.37 The rows depict the subdivision levels 1 through 3 of 7*. The red 
plots are the average y-value. The green plots are the average x,-values and the 
blue plots are the average x2-values. 
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Figure 3.38 Subdivision level 10 (the limit surface). 


Figure 3.39 The point cloud Dz colored red according to the average y-value, 
colored green according to the average x,-value, and then colored blue according to 
the average £2-value. 


to level 10, we get the idea of the limiting surface (Figure 3.38). When comparing 
the algorithm results Figure 3.38 and the actual point cloud Figure 3.39 we see that 
the algorithm result is very close to the actual point cloud with the exception of 
small jagged triangles along the edge of the circle. Most of the error in the average 


y—values occurs around the edges of the circle. 


Remarks 


In this section we consider a case that requires a more diverse selection of representa- 
tive measures to approximate the measure well. Let the measure p, of the reference 


triangle A*, be equivalent to the Lebesgue measure over the triangle 


A = ((—5/12, 0), (—10/51, 27/68), (1/3, 0). 
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Figure 3.40 The image on the left depicts the triangle with Lebesgue measure over 
A. The image on the right depicts the subdomain representative measure found by 
the algorithm. This poor representation leads to visible errors as we refine. 


Note that A is not contained in D and dp = y,4dx. Then the distribution-dependent 


quantities are approximately 


Moo(A*, p) = 14890, Myo(A*, p) = 0.013868, and Mo,(A*, p) = .0049786. 


Thus q(p) = (—0.093137, 0.033437) lies in the region T3 on the wv-plane, which corre- 
sponds to a Type 3 complementary measure. Following the calculations described in 
Section 3.3 for the computation of the complementary measure parameters, we have 
that b = .95463 and c = —.91433. Figure 3.40 displays both the measure p and the 
representative measure. The measure p is not approximated well by the subdomain 
representative measure. Furthermore, this (or similar) issue arises with any charac- 
teristic measure significantly different from our triangle measures and complementary 
measures. We are merely laying the foundation for work in this area and expected 
this type of error. Future work will explore higher order schemes which contain more 
diverse representative measures. 

The initial level of subdivision greatly impacts the limiting surface of the sub- 
division. If the measure of the initial level can be represented well by the triangle 
complementary measures in Section 3.3, then the limiting surface of the algorithm 
will also be good. On the other hand if the measure is poorly approximated initially, 


then that error will propagate through the refinement. Consider the uniform point 
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Figure 3.41 The first plot is the average y-values of the rectangular point cloud 
with 7*. The plot in the middle is the algorithm result after 5 additional levels of 
subdivision. The third plot is after 15 levels of subdivision. Notice the error around 
the corners of the rectangle. 


cloud on the rectangle 
R= ((-3, —A4), =3; 4), (3, 4), (3, —4)) 


with y = 1+2,—2. Then 4 triangles in 7* encounter the same corner issue discussed 
above, and that error is significant 5 levels of subdivision later and continues to the 
limit surface (Figure 3.41). Future work will explore higher order schemes which 


contain more diverse representative measures. 
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CHAPTER 4 


CONCLUSION 


In this dissertation we address the issues of irregularly distributed point clouds from 
two different perspectives. In Chapter 2 we address several critical issues arising in 
practical implementations of processing real point cloud data that exhibits irregu- 
larities. We examine some ideas known to work well for regular point clouds and 
modify them to extend their applicability for realistic raw data with a semi-regular 
underlying topology. We develop practical algorithms that realize certain ideas based 
on the theory and methodology in Learning Theory using adaptive partitioning from 
[4] and [2]. In particular, we address the issues of instability for higher order poly- 
nomial approximation described in [2] and suggest practical algorithms using trun- 
cated SVD to ensure that our approximation does not encounter this problem and 
to add to the stability. We employ our algorithm in the analysis and processing of 
LiDAR point clouds. In this realization, special criteria are developed for arriving 
at a nearly optimal adaptive partition for representation of the surface of interest; 
different adjustments to the results are applied to enhance the results and increase 
stability, including clustering to separate different sensed objects in the calculation 
of the approximation. A new nonlinear encoding process that progressively describes 
the surface approximations is designed for giving a realistic representation, especially 
for very large compression rates. 

Chapter 3 develops the theory and special algorithms targeted at representing 
curves and surfaces with gaps. These are particularly important for the representa- 


tion of the approximation results of the types of point cloud considered in Chapter 2. 


119 


While the standard subdivision schemes are developed to apply for regular distribu- 
tions, we suggest a new approach that gives the opportunity to extend this important 
visualization tool to new practically important areas. The algorithms analyze the ag- 
gregate quantities from the coarser level, separating them into distribution-dependent 
and value-dependent, to calculate their counterparts at a finer level. We give a ground- 
work for further theoretical developments in this area and the exploration of new 


distribution-dependent subdivision schemes. 


4.1 FUTURE WORK 


Here we describe our considerations for future work in the area. 

Related to the research discussed in Chapter 2 we have the following topics: mod- 
ification of the spherical learning algorithm to process point clouds from moving 
sensors; extension of the algorithm beyond d = 2 by partitioning with general sim- 
plices; and development of more sophisticated compression techniques to improve our 
rates. 

The ideas in Chapter 3 can be extended so that we have reproduction of higher 
order moments. In doing so, we will be able to have perfect recovery of more diverse 
types of measures. Along the same lines, a more sophisticated blending scheme could 
be implemented. We have already started exploring the polynomial of degree three 


weight measures and the corresponding subdomain measure. 
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